Searching for an suitable in vivo model of chronic ischemic tissue
Statistical report
Authors’ affiliations:
Veronika Lovasova1,2, Robert Bem3, Jaroslav Chlupac1,5, Michal Dubsky3,4, Jitka Husakova3,4, Andrea Nemcova3, Ludek Voska5, Zuzana Simunkova6, Jan Mares7, Filip Tichanek7, Jiri Fronek1,4,8
1Transplant Surgery Department, Institute for Clinical and Experimental Medicine, Prague, Czech Republic. 2Second Faculty of Medicine, Charles University, Prague, Czech Republic. 3Diabetic Center, Institute for Clinical and Experimental Medicine, Prague, Czech Republic. 4First Faculty of Medicine, Charles University, Prague, Czech Republic. 5Clinical and Transplant Pathology Centre, Institute for Clinical and Experimental Medicine, Prague, Czech Republic. 6Experimental Medicine Centre, Institute for Clinical and Experimental Medicine, Prague, Czech Republic. 7Department of Data Science, Institute for Clinical and Experimental Medicine, IKEM, Prague, Czech Republic. 8Department of Anatomy, Second Faculty of Medicine, Charles University, Prague, Czech Republic.
Person taking responsibility for the data analysis and code: Filip Tichanek (f.tichanek@gmail.com)
This is statistical report of the study that has been under review in the [TO BE ADDED]
Original publication
This is statistical report of for the [TO BE ADDED] published in the [TO BE ADDED]
When using this code, cite original publication:
TO BE ADDED
Original GitHub repository: https://filip-tichanek.github.io/wound_pig/
1 Introduction
Diabetes mellitus (DM) is associated with serious complications, such as diabetic foot ulcers (DFUs), wounds characterized by impaired healing. Peripheral artery disease (PAD), a major complication of DM, exacerbates the risk and severity of DFUs due to impaired tissue perfusion and local hypoxia, often leading to amputation. To address the urgent need for new treatments, our study focused on developing an in vivo model of chronic ischemic tissue that closely mimics non-healing ulcers in diabetic patients.
We introduced two models in diabetic and non-diabetic pigs, evaluating tissue oxygenation (TcPO2), wound healing rate (wound area change), and angiogenesis (Vessel index). We sought for model with characteristics indicating signs of chronicity, i.e. low and non-improving TcPO2 levels, delayed or absent wound closure and persistently impaired vascularization of the wound.
1.1 Study design
This study aimed to investigate the effects of ischemia on wound healing in both diabetic and non-diabetic pigs. Ischemic conditions were induced in limbs and through the creation of dorsal wounds, organizing the animals into four main groups based on their diabetic status and the localization of the wound:
- Limb wounds in non-diabetic pigs (
Non-DM, limbs) - Flap wounds in non-diabetic pigs (
Non-DM, flaps) - Limb wounds in diabetic pigs (
Diabetic, limbs) - Flap wounds in diabetic pigs (
Diabetic, flaps)
The experiment’s core evaluations included measurement of tissue oxygen levels (TcPO2), analysis of wound healing reflected by changes in wound area (or healing), and histological examination for vascular changes through the vessel index. Except for the vessel index, these assessments were conducted at several time points post-intervention.
TcPO2 (mmHg) was measured longitudinally: on the day of wound formation (D0), a week post-formation (D7), two weeks post-formation (D14), and a month after wound formation (D28).
Wound area (mm^2) was measured longitudinally on D0, D7, and D14.
Vessel index measurement requires a biopsy and was thus measured once per wound on D7 or D14 or D28.
1.2 Statistical methods
Due to the small sample size in our study, classical parametric statistics in a frequentist framework, which rely on large-sample approximations, were unsuitable. Thus, we opted for a Bayesian framework, which provides clear interpretations regardless of sample size (McElreath 2020).
We utilized hierarchical Bayesian regression models with either a Student t-distribution (‘robust’ model) or Gamma likelihood (with a log-link) in R software (R Core Team 2023), version 4.3.1. These models were employed to analyze the impact of diabetes (absent [non-DM] or present [diabetic]), wound site (on [flaps] or [limbs]), and their interaction (fitted via a summary variable group) on specified outcomes or their change over time. The models were fitted using the ‘brms’ package (Bürkner 2017, 2018), which employs ‘Stan’ software for probabilistic computing (Stan Development Team 2021).
There were three outcomes to be analyzed:
(i) Local oxygen measured via transcutaneous oximetry (
TcPO2, mmHg)(ii)
Wound area(mm^2)(iii)
Vessel index, reflecting tissue vascularization
To address data dependency, both the TcPO2 and wound area models incorporated hierarchical random effects. The random-effect factor for the wound was nested within the subject’s random effect. Given the low sample size, we primarily utilized a random intercept only. However, due to the presence of individual-specific trends in wound area, we also applied a more complex model with both random intercept and slope as secondary analysis.
Since the vessel index was measured only once per wound, this particular model included a simple random intercept for the subject.
The choice of likelihood distribution was guided by prior exploratory data analysis (EDA). Since TcPO2 and vessel index exhibited a right-tailed distribution with larger errors for larger predictions, a Gamma distribution with a log-link was chosen. Conversely, wound area did not demonstrate such characteristics, prompting us to model it using robust regression (Student t-distribution). These models were subsequently checked using a posterior predictive check (PPC) (McElreath 2020; Gelman et al. 2013).
To ensure all models fit the data well and met distributional assumptions, we conducted the PPC for each model. Additionally, we examined Pareto k values to detect overly influential or outlying data points (Vehtari et al. 2020). To compare the predictive accuracy of different models, we employed leave-one-out cross-validation (LOO-CV) (Vehtari, Gelman, and Gabry 2016).
Conservative priors centered at zero effect were generally used for the group and group:time interaction effects. Priors were Gaussian, with \(mu = 0\) and \(sigma = 2\) (for robust model of wound area, \(sigma = 2*SD\) where \(SD\) being standard deviation of wound area) for the effect group. Time was scaled to the unit of weeks. For the time:group interaction, \(sigma = 1\) was used (1 standard deviation of the wound area for robust model respectively).
For the effect of time, we assumed that TcPO2 should increase over time, reflecting healing. Similarly, according to Roy et al. (Roy et al. 2008), the wound area should decrease by at least 20% per 14 days, i.e., approximately by 50 mm2 per week in the context of our data. This assumption was reflected in the prior distribution for the time effect, with:
\(normal (mu = 0.01, sigma = 2)\) for Gamma models
\(normal (mu = -50, sigma = 350)\) for the robust model (350 = 2 standard deviations of the
wound area)
To report the uncertainty of the effect estimation, we show 95% credible intervals and the probability of direction (PD), which can be interpreted as an index (0.5 to 1) representing the certainty that the effect goes in a specific direction (Makowski et al. 2019).
Data, statistical code, and a detailed analytical report are publicly available on https://github.com/filip-tichanek/pigWound from the day of publication.
2 Analysis
2.1 Data and packages
2.1.1 Packages
Open code
if (T) {rm(list = ls() )}
getwd()
## [1] "/home/ticf/GitRepo/ticf/190_LOVV_diabetes_wound"
if (T) {
suppressWarnings(suppressMessages({
library(tidyverse)
library(stringr)
library(emmeans)
library(gtsummary)
library(car)
library(sjPlot)
library(flextable)
library(openxlsx)
library(mgcv)
library(pROC)
library(cowplot)
library(boot)
library(glmnet)
library(brms)
library(projpred)
library(janitor)
library(arm)
library(corrplot)
library(bayesplot)
library(ggbeeswarm)
library(kableExtra)
# Functions clashes
select <- dplyr::select
rename <- dplyr::rename
mutate <- dplyr::mutate
recode <- dplyr::recode
summarize <- dplyr::summarize
count <- dplyr::count
# Simple functions
logit <-function(x){log(x/(1-x))}
inv.logit <- function(x){exp(x)/(1+exp(x))}
}))
}2.1.2 Functions
2.1.2.1 run_model
The function to (i) load OR (ii) run & save (if has not been done previously) results of any computation or complex table production
Open code
run_model <- function(expr, path, reuse = TRUE) {
# Initialize 'fit' to NULL to ensure it's always defined
fit <- NULL
# Construct the file path only if 'reuse' is TRUE
if (reuse) {
path <- paste0(path, ".Rds")
fit <- suppressWarnings(try(readRDS(path), silent = TRUE))
# Check if 'fit' is an error (file not found or couldn't be read)
if (inherits(fit, "try-error")) {
fit <- NULL
}
}
# If 'fit' is NULL (either because 'reuse' is FALSE,
# or the file couldn't be read), evaluate 'expr'
if (is.null(fit)) {
fit <- eval(substitute(expr))
# Save 'fit' only if 'reuse' is TRUE and a valid 'path' is provided
if (reuse && !is.null(path) && nzchar(path)) {
saveRDS(fit, file = path)
}
}
return(fit)
}2.1.2.2 p_dir
The function calculates proportion of posterior probability behind specified threshold or probability of direction (as defined in Makowski et al., 2019 (Makowski et al. 2019)). data: data frame containing posterior draws in columns. dir: probability under/above threshold (> or <) or probability of direction (max). threshold: often null effect
Open code
p_dir <- function(data, dir, tres){
if(dir == 'max'){
return(
max(length(data[data>tres]),length(data[data<tres]))/length(data)
)
} else if(dir == '<'){
return(
length(data[data<tres])/length(data)
)
} else if(dir == '>'){
return(
length(data[data>tres])/length(data)
)
} else {print('ERROR')}
}2.1.2.3 cur
The function creates polygon or curve given the data frame of posterior samples (post), specific column given (var), position on X axis (con), filling color (cole) and color of border of the polygon (borde). The last two arguments have defaults.
Open code
cur <- function(post, var, con, cole, borde){
de <- density(post[[var]])[c('x','y')] %>% data.frame()
if(missing(cole)) {cole = rgb(1, .5, .1, alpha=0.4)}
if(missing(borde)) {borde = cole}
polygon(x = de$x, y = de$y + con, col = cole, border = borde)
}2.1.2.4 ci
The function draws credible intervals (95% in wide line and 99% as slimmer line), given the data frame of posterior samples (post), specified column (var), position on X axis (con), and color cole. It also automatically generates white line on the same Y-axis as is the placement of the CI
Open code
ci <- function(post, var, con, cole, xsc){
if(missing(cole)) {cole = rgb(.9, .4, 0)}
qs = quantile(post[[var]], probs = c(0.005, 0.995, 0.025, 0.975, 0.5))
lines(c(min(xsc), max(xsc) ), c(con, con), col = 'white')
lines(c(qs[1], qs[2]), c(con, con), col = cole, lwd=1.3, lend = 1)
lines(c(qs[3], qs[4]), c(con, con), col = cole, lwd=3.5, lend=1)
points(qs[5], con, pch='I', cex = 1)
}Setting directories
Open code
## directories
if (file.exists('data') == FALSE) {
DATA_DIR <- "data"
dir.create(DATA_DIR, showWarnings = FALSE, recursive = TRUE)
}
if (file.exists('output') == FALSE) {
OUTPUT_DIR <- "gitignore/output"
dir.create(OUTPUT_DIR, showWarnings = FALSE, recursive = TRUE)
}
if (file.exists('output/brm') == FALSE) {
BRM_DIR <- "gitignore/output/brm"
dir.create(BRM_DIR, showWarnings = FALSE, recursive = TRUE)
}
if (file.exists('data/db_history') == FALSE) {
DB_HISTORY_DIR <- "data/db_history"
dir.create(DB_HISTORY_DIR, showWarnings = FALSE, recursive = TRUE)
}2.1.3 Import data
Open code
dat <- read.xlsx("data/data_Lovasova_tidy_deleted.xlsx")
dat$Subject <- as.factor(dat$Subject)
dat$Site <- as.factor(dat$Site)
dat$Type <- as.factor(dat$Type)
dat$Site_nest <- interaction(dat$Subject, dat$Site)
dat$DMcat <- ifelse(dat$DM == 1, 'DM1', 'DM0')
dat$group <- interaction(dat$DMcat, dat$Type)
dat$group <- factor(dat$group, levels = c('DM0.flaps', 'DM0.limbs', 'DM1.flaps', 'DM1.limbs'))
dat$group2 <- as.factor(ifelse(dat$Control == 1, 'in_Control', as.character(dat$group) ) )
summary(dat)
## Subject Site DM Type Control Vessel_7
## 1 : 5 1:12 Min. :0.0 flaps:30 Min. :0.0 Min. : 38.00
## 2 : 5 2:12 1st Qu.:0.0 limbs:30 1st Qu.:0.0 1st Qu.: 65.00
## 3 : 5 3:12 Median :0.5 Median :0.0 Median : 68.00
## 4 : 5 4:12 Mean :0.5 Mean :0.4 Mean : 81.35
## 5 : 5 5:12 3rd Qu.:1.0 3rd Qu.:1.0 3rd Qu.:108.00
## 6 : 5 Max. :1.0 Max. :1.0 Max. :132.00
## (Other):30 NA's :43
## Vessel_14 Vessel_28 Healing_0 Healing_7
## Min. : 39.00 Min. : 37.00 Min. :119.0 Min. :131.0
## 1st Qu.: 54.50 1st Qu.: 54.50 1st Qu.:372.2 1st Qu.:314.2
## Median : 74.00 Median : 64.00 Median :443.0 Median :418.5
## Mean : 72.68 Mean : 66.75 Mean :456.4 Mean :411.6
## 3rd Qu.: 84.00 3rd Qu.: 74.50 3rd Qu.:539.0 3rd Qu.:498.5
## Max. :128.00 Max. :110.00 Max. :844.0 Max. :792.0
## NA's :41 NA's :40 NA's :4
## Healing_14 TcPO2_0 TcPO2_7 TcPO2_14
## Min. : 27.0 Min. : 0.00 Min. : 1.00 Min. : 1.00
## 1st Qu.:102.5 1st Qu.: 4.50 1st Qu.: 4.75 1st Qu.: 5.75
## Median :184.5 Median :15.00 Median :16.50 Median :20.00
## Mean :237.4 Mean :25.07 Mean :23.65 Mean :22.83
## 3rd Qu.:335.0 3rd Qu.:44.50 3rd Qu.:44.00 3rd Qu.:35.00
## Max. :883.0 Max. :81.00 Max. :69.00 Max. :61.00
## NA's :1
## TcPO2_28 Site_nest DMcat group group2
## Min. : 1.00 1.1 : 1 Length:60 DM0.flaps:15 DM0.flaps :12
## 1st Qu.:15.00 2.1 : 1 Class :character DM0.limbs:15 DM0.limbs : 6
## Median :29.00 3.1 : 1 Mode :character DM1.flaps:15 DM1.flaps :12
## Mean :28.04 4.1 : 1 DM1.limbs:15 DM1.limbs : 6
## 3rd Qu.:41.00 5.1 : 1 in_Control:24
## Max. :60.00 6.1 : 1
## NA's :5 (Other):542.1.4 Wrangle data
In the next step, we will modify data as needed for exploration and modelling. Note that we also create a variable TcPO0_nz where we replace zero values with 0.5. This is to address zeros that arise from detection limits, not actual absence, “true zero”. This enables fitting Gamma models that do not allow occurrence of zeros.
2.1.4.1 TcPO2
Open code
tcPO_long_all <- dat %>% pivot_longer(
cols = starts_with('Tc'),
names_to = 'time_days',
values_to = 'TcPO2'
) %>% select(
Subject, Site, Site_nest, DM, Type, group, Control, time_days, TcPO2, group2
) %>% filter(!is.na(TcPO2)) %>% mutate(
time_days = as.numeric(
gsub('TcPO2_', '', time_days)
),
DM = if_else(DM == 0, -0.5, 0.5),
limbs = if_else(Type == 'limbs', 0.5, -0.5),
TcPO2_nz = if_else(TcPO2 == 0, 0.5, TcPO2),
time_cen = (time_days - 14) / 7
)
tcPO_long <- tcPO_long_all %>% filter(Control == 0)
summary(tcPO_long)
## Subject Site Site_nest DM Type group
## 2 :16 1:47 2.1 : 4 Min. :-0.50000 flaps:91 DM0.flaps:44
## 3 :16 2:47 3.1 : 4 1st Qu.:-0.50000 limbs:48 DM0.limbs:24
## 7 :16 3:23 4.1 : 4 Median : 0.50000 DM1.flaps:47
## 8 :16 4:22 5.1 : 4 Mean : 0.01079 DM1.limbs:24
## 9 :15 5: 0 6.1 : 4 3rd Qu.: 0.50000
## 1 :12 7.1 : 4 Max. : 0.50000
## (Other):48 (Other):115
## Control time_days TcPO2 group2 limbs
## Min. :0 Min. : 0.00 Min. : 0.00 DM0.flaps :44 Min. :-0.5000
## 1st Qu.:0 1st Qu.: 3.50 1st Qu.: 3.00 DM0.limbs :24 1st Qu.:-0.5000
## Median :0 Median : 7.00 Median : 8.00 DM1.flaps :47 Median :-0.5000
## Mean :0 Mean :11.88 Mean :13.22 DM1.limbs :24 Mean :-0.1547
## 3rd Qu.:0 3rd Qu.:14.00 3rd Qu.:18.00 in_Control: 0 3rd Qu.: 0.5000
## Max. :0 Max. :28.00 Max. :60.00 Max. : 0.5000
##
## TcPO2_nz time_cen
## Min. : 0.50 Min. :-2.0000
## 1st Qu.: 3.00 1st Qu.:-1.5000
## Median : 8.00 Median :-1.0000
## Mean :13.23 Mean :-0.3022
## 3rd Qu.:18.00 3rd Qu.: 0.0000
## Max. :60.00 Max. : 2.0000
## 2.1.4.2 Wound area
Open code
healing_long_all <- dat %>% pivot_longer(
cols = starts_with('Heal'),
names_to = 'time_days',
values_to = 'healing'
) %>% select(
Subject, Site, Site_nest, DM, Type, group, Control, time_days, healing, group2
) %>% filter(!is.na(healing)) %>% mutate(
time_days = as.numeric(
gsub('Healing_', '', time_days)
),
DM = if_else(DM == 0, -0.5, 0.5),
limbs = if_else(Type == 'limbs', 0.5, -0.5),
time_cen = (time_days - 7)/7
)
healing_long <- healing_long_all %>% filter(Control == 0)
summary(healing_long)
## Subject Site Site_nest DM Type group
## 1 :12 1:36 1.1 : 3 Min. :-0.50000 flaps:69 DM0.flaps:36
## 2 :12 2:35 2.1 : 3 1st Qu.:-0.50000 limbs:36 DM0.limbs:18
## 3 :12 3:17 3.1 : 3 Median :-0.50000 DM1.flaps:33
## 7 :12 4:17 4.1 : 3 Mean :-0.01429 DM1.limbs:18
## 9 :12 5: 0 5.1 : 3 3rd Qu.: 0.50000
## 8 : 9 6.1 : 3 Max. : 0.50000
## (Other):36 (Other):87
## Control time_days healing group2 limbs
## Min. :0 Min. : 0 Min. : 39.0 DM0.flaps :36 Min. :-0.5000
## 1st Qu.:0 1st Qu.: 0 1st Qu.:245.0 DM0.limbs :18 1st Qu.:-0.5000
## Median :0 Median : 7 Median :395.0 DM1.flaps :33 Median :-0.5000
## Mean :0 Mean : 7 Mean :374.7 DM1.limbs :18 Mean :-0.1571
## 3rd Qu.:0 3rd Qu.:14 3rd Qu.:493.0 in_Control: 0 3rd Qu.: 0.5000
## Max. :0 Max. :14 Max. :883.0 Max. : 0.5000
##
## time_cen
## Min. :-1
## 1st Qu.:-1
## Median : 0
## Mean : 0
## 3rd Qu.: 1
## Max. : 1
## 2.1.4.3 Vessel index
Open code
vessel_long <- dat %>% pivot_longer(
cols = starts_with('Vess'),
names_to = 'time_days',
values_to = 'vessel_index'
) %>% select(
Subject, Site, Site_nest, DM, Type, group, Control, time_days, vessel_index, group2
) %>% filter(!is.na(vessel_index),
Control == 0) %>% mutate(
time_days = as.numeric(
gsub('Vessel_', '', time_days)
),
DM = if_else(DM == 0, -0.5, 0.5),
limbs = if_else(Type == 'limbs', 0.5, -0.5),
time_cen = (time_days - 17.5) / 7
)
summary(vessel_long)
## Subject Site Site_nest DM Type group
## 2 : 4 1:12 1.1 : 1 Min. :-0.50000 flaps:22 DM0.flaps:10
## 3 : 4 2:11 2.1 : 1 1st Qu.:-0.50000 limbs:12 DM0.limbs: 6
## 7 : 4 3: 5 3.1 : 1 Median : 0.50000 DM1.flaps:12
## 8 : 4 4: 6 4.1 : 1 Mean : 0.02941 DM1.limbs: 6
## 9 : 4 5: 0 5.1 : 1 3rd Qu.: 0.50000
## 1 : 2 6.1 : 1 Max. : 0.50000
## (Other):12 (Other):28
## Control time_days vessel_index group2 limbs
## Min. :0 Min. : 7.00 Min. : 43.00 DM0.flaps :10 Min. :-0.5000
## 1st Qu.:0 1st Qu.: 7.00 1st Qu.: 58.00 DM0.limbs : 6 1st Qu.:-0.5000
## Median :0 Median :14.00 Median : 71.00 DM1.flaps :12 Median :-0.5000
## Mean :0 Mean :16.88 Mean : 74.59 DM1.limbs : 6 Mean :-0.1471
## 3rd Qu.:0 3rd Qu.:28.00 3rd Qu.: 83.50 in_Control: 0 3rd Qu.: 0.5000
## Max. :0 Max. :28.00 Max. :132.00 Max. : 0.5000
##
## time_cen
## Min. :-1.50000
## 1st Qu.:-1.50000
## Median :-0.50000
## Mean :-0.08824
## 3rd Qu.: 1.50000
## Max. : 1.50000
## 2.2 Data exploration
2.2.1 Data distributions of outcomes, SupFig1
Open code
par(mfrow = c(1, 3))
hist((tcPO_long$TcPO2), 14, main = 'TcPO2', xlab='')
hist((healing_long$healing), 14, main = 'wound area', xlab='')
hist((vessel_long$vessel_index), 14, main = 'vessel index', xlab='')
2.2.2 TcPO2
2.2.2.1 All subjects per plot
Each plot include individuals from the same group (defined by presence of diabetes and location of wound)
Open code
cole <- c("#FF0000", "#80AF00", "#8000FF",
"#FF8000", "#01AF40", "#0000FF",
"#FF0080", "#90BF51", "#0001FF",
"#FF00FF", "#00CF00", "#0010FF")
fig_01a <- tcPO_long %>%
ggplot(aes(x= time_days, y=TcPO2, by = Site_nest, col = Subject)) +
geom_line(alpha = 0.8) +
facet_wrap(~group, ncol=4, labeller = labeller(group = c("DM0.flaps" = "Non-DM, flaps",
"DM0.limbs" = "Non-DM, limbs",
"DM1.flaps" = "Diabetic, flaps",
"DM1.limbs" = "Diabetic, limbs") ) ) +
scale_color_manual(values = cole) +
labs(x = "Days of experiment", y = 'TcPO2 (mmHg)') +
scale_x_continuous(breaks=c(0, 7, 14, 28))+
guides(color = guide_legend(override.aes = list(size = 3, linewidth=1.1), ncol=2))
fig_01a
2.2.2.2 One subject per plot
Each individual has a separate plots. Groups are distinguishable by color.
Open code
cole <- c("#01AF40", "#0000FF","#FF0000",
"#FF00FF", 'grey60')
fig_01b <- tcPO_long_all %>%
ggplot(aes(x = time_days, y = TcPO2, by = Site_nest, col = group2)) +
geom_line(alpha = 0.9) +
geom_point(alpha = 0.9) +
facet_wrap(~ Subject, ncol = 4, dir = 'v',
labeller = labeller(Subject = function(x) paste0("Subject ", x))) +
scale_color_manual(values=cole,
name="Group",
breaks=c("DM0.flaps", "DM0.limbs", "DM1.flaps",
"DM1.limbs", "in_Control"),
labels=c("Non-DM, flaps", "Non-DM, limbs", "Diabetic, flaps", "Diabetic, limbs",
"Control")) +
labs(x = "Days of experiment", y = 'TcPO2 (mmHg)') +
scale_x_continuous(breaks=c(0, 7, 14, 28))+
theme(legend.position="bottom",
legend.key.size = unit(0.4, 'cm'),
legend.title = element_text(size=11),
legend.text = element_text(size=10))+
guides(color = guide_legend(override.aes = list(size = 2, linewidth=1), ncol=4))
fig_01b
2.2.3 Wound area
2.2.3.1 All subjects per plot
Each plot include individuals from the same group (defined by presence of diabetes and location of wound)
Open code
cole <- c("#FF0000", "#80AF00", "#8000FF",
"#FF8000", "#01AF40", "#0000FF",
"#FF0080", "#90BF51", "#0001FF",
"#FF00FF", "#00CF00", "#0010FF")
fig_02a <- healing_long %>%
ggplot(aes(x= time_days, y=healing, by = Site_nest, col = Subject)) +
geom_line(alpha = 0.8) +
facet_wrap(~group, ncol=4, labeller = labeller(group = c("DM0.flaps" = "Non-DM, flaps",
"DM0.limbs" = "Non-DM, limbs",
"DM1.flaps" = "Diabetic, flaps",
"DM1.limbs" = "Diabetic, limbs") ) ) +
scale_color_manual(values = cole) +
labs(x = "Days of experiment", y = expression(Wound~area~(mm^2))) +
scale_x_continuous(breaks=c(0, 7, 14))+
guides(color = guide_legend(override.aes = list(size = 3, linewidth=1.1), ncol=2))
fig_02a
2.2.3.2 One subject per plot
Each individual has a separate plots. Groups are distinguishable by color.
Open code
cole <- c("#01AF40", "#0000FF","#FF0000",
"#FF00FF", 'grey60')
fig_02b <- healing_long_all %>%
ggplot(aes(x = time_days, y = healing, by = Site_nest, col = group2)) +
geom_line(alpha = 0.9) +
geom_point(alpha = 0.9) +
facet_wrap(~ Subject, ncol = 4, dir = 'v',
labeller = labeller(Subject = function(x) paste0("Subject ", x))) +
scale_color_manual(values=cole,
name="Group",
breaks=c("DM0.flaps", "DM0.limbs", "DM1.flaps",
"DM1.limbs", "in_Control"),
labels=c("Non-DM, flaps", "Non-DM, limbs", "Diabetic, flaps", "Diabetic, limbs",
"Control")) +
labs(x = "Days of experiment", y = expression(Wound~area~(mm^2))) +
scale_x_continuous(breaks=c(0, 7, 14))+
theme(legend.position="bottom",
legend.key.size = unit(0.4, 'cm'),
legend.title = element_text(size=11),
legend.text = element_text(size=10))+
guides(color = guide_legend(override.aes = list(size = 2, linewidth=1), ncol=4))
fig_02b
2.2.4 Vessel index
Note that there is only a single measurement per wound for the vessel index
Each plot include individuals from the same group (defined by presence of diabetes and location of wound). Individuals are distinguishable by color.
Open code
cole <- c("#FF0000", "#80AF00", "#8000FF",
"#FF8000", "#01AF40", "#0000FF",
"#FF0080", "#90BF51", "#0001FF",
"#FF00FF", "#00CF00", "#0010FF")
fig_03a <- vessel_long %>%
ggplot(aes(x= time_days, y=vessel_index, col = Subject)) +
geom_beeswarm(alpha = 0.8, size = 3) +
scale_color_manual(values = cole) +
stat_smooth(data=subset(vessel_long, group=="DM0.flaps"),
aes(x=time_days, y=vessel_index),
method = 'lm', se = F, color='black', size = 0.8) +
stat_smooth(data=subset(vessel_long, group=="DM0.limbs"),
aes(x=time_days, y=vessel_index),
method = 'lm', se = F, color='black', size = 0.8) +
stat_smooth(data=subset(vessel_long, group=="DM1.flaps"),
aes(x=time_days, y=vessel_index),
method = 'lm', se = F, color='black', size = 0.8) +
stat_smooth(data=subset(vessel_long, group=="DM1.limbs"),
aes(x=time_days, y=vessel_index),
method = 'lm', se = F, color='black', size = 0.8) +
facet_grid(~group, labeller = labeller(group = c(
"DM0.flaps" = "Non-DM, flaps",
"DM0.limbs" = "Non-DM, limbs",
"DM1.flaps" = "Diabetic, flaps",
"DM1.limbs" = "Diabetic, limbs") ) ) +
labs(x = "Days of experiment", y = 'Vessel index') +
theme(legend.position="bottom",
legend.key.size = unit(0.4, 'cm'),
legend.title = element_text(size=14),
legend.text = element_text(size=12)) +
guides(color = guide_legend(override.aes = list(size = 4.6, linewidth=1), nrow=2))+
scale_x_continuous(breaks=c(7, 14, 28))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
fig_03a
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
2.3 Models
2.3.1 TcPO2 model
Because TcPO2 is only-positive variable with heavy right tail (Supplementary Figure 1), we will primarily model it via hierarchical generalized linear model with Gamma distribution and log-link.
2.3.1.1 Priors
We did not find any data about development of TcPO2 over time in similar pig model. However, we still expect the TcPO2 will be at least slightly improving over time. Thus, the prior will reflect the expectation of increasing the TcPO2 1.1 times per a week, thereby setting the center of prior for the time effect slightly
Open code
log(1.1)
## [1] 0.09531018Setting the prior distribution
Open code
prior1 <- prior(normal(0, 2), class= b, coef = 'groupDM0.limbs') +
prior(normal(0, 2), class= b, coef = 'groupDM1.flaps')+
prior(normal(0, 2), class= b, coef = 'groupDM1.limbs')+
prior(normal(0.01, 2), class= b, coef = 'time_cen')+
prior(normal(0, 1), class= b, coef = 'time_cen:groupDM0.limbs')+
prior(normal(0, 1), class= b, coef = 'time_cen:groupDM1.flaps')+
prior(normal(0, 1), class= b, coef = 'time_cen:groupDM1.limbs')2.3.1.2 Model fit
Next, we will fit the model and also evaluate model convergence and effective size of posterior samples (ESS). We will use following parameters:
Rhatfrom model summary: 1 indicates convergenceBulk_ESSandTail_ESSfrom summary output: should be > 1000-2000, ideally more for effects of interest
Open code
tcpo_model_int <- run_model(
brm(TcPO2_nz ~ time_cen + group + time_cen:group + (1|Subject/Site_nest),
data = tcPO_long,
prior = prior1,
chains = 2, iter = 8000, warmup = 2000, seed = 17,
control = list(adapt_delta = 0.99),
cores = 2,
family = Gamma(link='log')),
'gitignore/output/brm/tcpo_model_int', reuse = TRUE)
summary(tcpo_model_int)
## Family: gamma
## Links: mu = log; shape = identity
## Formula: TcPO2_nz ~ time_cen + group + time_cen:group + (1 | Subject/Site_nest)
## Data: tcPO_long (Number of observations: 139)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Subject (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.43 0.20 0.09 0.90 1.00 2653 2668
##
## ~Subject:Site_nest (Number of levels: 36)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.18 0.12 0.01 0.45 1.00 3583 5784
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 1.98 0.30 1.42 2.63 1.00 4793
## time_cen 0.30 0.11 0.09 0.51 1.00 4550
## groupDM0.limbs 0.98 0.44 0.05 1.82 1.00 5985
## groupDM1.flaps 0.61 0.42 -0.29 1.41 1.00 5112
## groupDM1.limbs 0.94 0.44 0.01 1.77 1.00 5275
## time_cen:groupDM0.limbs -0.06 0.17 -0.39 0.27 1.00 6986
## time_cen:groupDM1.flaps -0.27 0.14 -0.54 0.01 1.00 5811
## time_cen:groupDM1.limbs 0.08 0.16 -0.24 0.40 1.00 7151
## Tail_ESS
## Intercept 5020
## time_cen 6892
## groupDM0.limbs 6459
## groupDM1.flaps 5308
## groupDM1.limbs 6095
## time_cen:groupDM0.limbs 9010
## time_cen:groupDM1.flaps 8814
## time_cen:groupDM1.limbs 8922
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 1.26 0.15 0.99 1.56 1.00 13035 9190
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Convergence of model is fine, ESSs are sufficient.
Lets look at estimated trend in TcPO2 over time and across groups
Open code
conditional_effects(tcpo_model_int, effects='time_cen:group')
Now look below to predicted TcPO2 in the middle day, separately for each group
Open code
conditional_effects(tcpo_model_int, effects='group')
Based on both plots, it seems that wounds on flaps may restore TcPO2 more slowly compared to limbs
2.3.1.3 Model diagnostics
These plots show how well model reproduce (predict, or retrodict) data
We will show following plot:
- data density plot: the figure shows density of the outcome variable (thick line) and densities of artificial data from simulations based on Bayesian model estimates (slimmer lines). The thick line should not deviate much from the slimmer lines (in that case, model does not recapitulate data well). We will show specifically grouped data density which is the same but for specific group. This is also useful for evaluation potential group-specific shapes of data distribution (e.g. different spread of outcome)
Open code
pp_check(tcpo_model_int, type='dens_overlay_grouped',ndraws = 50,group='group') + xlim(0, 100)
It looks very fine here as thick line closely follows the slim lines.
- scatter plot shows average predictions (x-axis) vs. data. For Gaussian and robust (t-distribution) models, increasing X-axis should not noticeably change the spread of data points
Open code
pp_check(tcpo_model_int, type='scatter_avg')
As expected, spread increases with X axis. This would be problem if we had Gaussian model. As we have Gamma model, this behavior is expected,
- max - min plot shows predicted vs. actual minimal and maximal value of outcome. Dark blue point should be among the light-blue points
Open code
set.seed(16)
pp_check(tcpo_model_int,type="stat_2d", stat = c("max", "min"),ndraws=200)
This is not good here. Our minimal value is relatively large (this is likely product of our addition of to zero) and our maximal value is surprisingly small. TcPO2 thus seems to have upper bound, limiting its further increase.
- mean - standard deviation plot shows predicted vs. actual mean and SD of the outcome
Open code
set.seed(16)
pp_check(tcpo_model_int,type="stat_2d", stat = c("mean", "sd"),ndraws=200)
This seems relatively fine: the model expects the mean and standard deviation of the TcPO2 not far away from the factual values.
- LOO plots answers the question ‘are there too influential observations?’ and show calibration plot respectively. All pareto estimates should have k < 0.7 (indicating no overly influential observation) and points in the second plot should be close the line, indicating a good calibration.
Open code
plot(loo(tcpo_model_int))
Open code
pp_check(tcpo_model_int, type = "loo_pit")
Open code
loo(tcpo_model_int)
##
## Computed from 12000 by 139 log-likelihood matrix
##
## Estimate SE
## elpd_loo -484.9 14.4
## p_loo 18.3 3.0
## looic 969.8 28.8
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 129 92.8% 2168
## (0.5, 0.7] (ok) 8 5.8% 561
## (0.7, 1] (bad) 2 1.4% 271
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.There are few influential observations as would be expected given small sample size. We have to keep in mind that model may be overconfident and the estimated uncertainty underestimated.
2.3.2 Wound area models
The wound area, unlike the other two outcomes, displays a symmetrical distribution, suggesting its suitability for Gaussian approximation. However, due to the small sample size, a robust regression with a Student’s t-distribution is recommended to mitigate the impact of outliers.
Figure of changes in wound area over time (the section 2.2.3. above) shows considerable individual variation even within a single group, indicating the need for a more sophisticated model. For our secondary analysis, we will use a model with random both intercepts and slopes to accurately account for individual and temporal variations.
2.3.2.1 Prior specification
According to Roy et al., the wound area should decrease by at least 20% per 14 days, i.e. approximately by 50 units per week in the context of our data. This will be addressed during the prior setting.
Data summary for prior specification
Open code
sd(healing_long$healing) * c(2, 1)
## [1] 350.1426 175.0713Prior specification
Open code
prior2 <- prior(normal(0, 350), class= b, coef = 'groupDM0.limbs') +
prior(normal(0, 350), class= b, coef = 'groupDM1.flaps')+
prior(normal(0, 350), class= b, coef = 'groupDM1.limbs')+
prior(normal(-50, 350), class= b, coef = 'time_cen')+
prior(normal(0, 175), class= b, coef = 'time_cen:groupDM0.limbs')+
prior(normal(0, 175), class= b, coef = 'time_cen:groupDM1.flaps')+
prior(normal(0, 175), class= b, coef = 'time_cen:groupDM1.limbs')2.3.2.2 Random intercept model - Gaussian
Open code
healing_model_int <- run_model(
brm(healing ~ time_cen + group + time_cen:group + (1|Subject/Site_nest),
data = healing_long,
prior = prior2,
chains = 2, iter = 8000, warmup = 2000, seed = 17,
control = list(adapt_delta = 0.99),
cores = 2,
family = gaussian() ,
save_pars = save_pars(all = TRUE)),
'gitignore/output/brm/healing_model_int', reuse = TRUE)
summary(healing_model_int)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: healing ~ time_cen + group + time_cen:group + (1 | Subject/Site_nest)
## Data: healing_long (Number of observations: 105)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Subject (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 34.69 25.11 1.47 95.01 1.00 3075 4411
##
## ~Subject:Site_nest (Number of levels: 36)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 41.49 22.68 2.61 86.31 1.00 2627 4351
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 417.14 34.20 347.22 483.87 1.00 7331
## time_cen -64.27 25.10 -113.22 -14.65 1.00 8317
## groupDM0.limbs -151.77 54.62 -260.28 -43.98 1.00 8180
## groupDM1.flaps -71.09 49.26 -167.46 28.84 1.00 8096
## groupDM1.limbs 32.33 54.21 -74.87 140.02 1.00 7873
## time_cen:groupDM0.limbs -4.84 43.50 -90.02 80.63 1.00 10850
## time_cen:groupDM1.flaps -75.12 35.31 -145.53 -6.21 1.00 9525
## time_cen:groupDM1.limbs -103.14 43.41 -187.87 -16.43 1.00 10594
## Tail_ESS
## Intercept 6135
## time_cen 9155
## groupDM0.limbs 7652
## groupDM1.flaps 7317
## groupDM1.limbs 6860
## time_cen:groupDM0.limbs 9049
## time_cen:groupDM1.flaps 9589
## time_cen:groupDM1.limbs 8845
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 128.46 10.81 108.73 151.40 1.00 6797 8668
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(healing_model_int, effects='time_cen:group')
Open code
conditional_effects(healing_model_int, effects='group')
Model converged well and shows sufficient ESS. The model suggests that wounds of diabetic pigs may have, surprisingly, faster wound closure. Very similar result are expected also from robust model
2.3.2.3 Random intercept model - robust
Open code
healing_model_int_robust <- run_model(
brm(healing ~ time_cen + group + time_cen:group + (1|Subject/Site_nest),
data = healing_long,
prior = prior2,
chains = 2, iter = 8000, warmup = 2000, seed = 17,
control = list(adapt_delta = 0.99),
cores = 2,
family = student() ,
save_pars = save_pars(all = TRUE)),
'gitignore/output/brm/healing_model_int_robust', reuse = TRUE)
summary(healing_model_int)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: healing ~ time_cen + group + time_cen:group + (1 | Subject/Site_nest)
## Data: healing_long (Number of observations: 105)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Subject (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 34.69 25.11 1.47 95.01 1.00 3075 4411
##
## ~Subject:Site_nest (Number of levels: 36)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 41.49 22.68 2.61 86.31 1.00 2627 4351
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 417.14 34.20 347.22 483.87 1.00 7331
## time_cen -64.27 25.10 -113.22 -14.65 1.00 8317
## groupDM0.limbs -151.77 54.62 -260.28 -43.98 1.00 8180
## groupDM1.flaps -71.09 49.26 -167.46 28.84 1.00 8096
## groupDM1.limbs 32.33 54.21 -74.87 140.02 1.00 7873
## time_cen:groupDM0.limbs -4.84 43.50 -90.02 80.63 1.00 10850
## time_cen:groupDM1.flaps -75.12 35.31 -145.53 -6.21 1.00 9525
## time_cen:groupDM1.limbs -103.14 43.41 -187.87 -16.43 1.00 10594
## Tail_ESS
## Intercept 6135
## time_cen 9155
## groupDM0.limbs 7652
## groupDM1.flaps 7317
## groupDM1.limbs 6860
## time_cen:groupDM0.limbs 9049
## time_cen:groupDM1.flaps 9589
## time_cen:groupDM1.limbs 8845
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 128.46 10.81 108.73 151.40 1.00 6797 8668
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(healing_model_int_robust, effects='time_cen:group')
Open code
conditional_effects(healing_model_int_robust, effects='group')
Again, model converged well and shows sufficient ESS. Results are very similar to Gaussian model as expected. Will more complex and presumably more conservative random-slope model replicate these findings?
2.3.2.4 Random slope model
Open code
healing_model_int_rslope <- run_model(
brm(healing ~ time_cen + group + time_cen:group + (time_cen|Subject/Site_nest),
data = healing_long,
prior = prior2,
chains = 2, iter = 8000, warmup = 2000, seed = 17,
control = list(adapt_delta = 0.99),
cores = 2,
family = student() ,
save_pars = save_pars(all = TRUE)),
'gitignore/output/brm/healing_model_int_rslope', reuse = TRUE)
summary(healing_model_int)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: healing ~ time_cen + group + time_cen:group + (1 | Subject/Site_nest)
## Data: healing_long (Number of observations: 105)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Subject (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 34.69 25.11 1.47 95.01 1.00 3075 4411
##
## ~Subject:Site_nest (Number of levels: 36)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 41.49 22.68 2.61 86.31 1.00 2627 4351
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 417.14 34.20 347.22 483.87 1.00 7331
## time_cen -64.27 25.10 -113.22 -14.65 1.00 8317
## groupDM0.limbs -151.77 54.62 -260.28 -43.98 1.00 8180
## groupDM1.flaps -71.09 49.26 -167.46 28.84 1.00 8096
## groupDM1.limbs 32.33 54.21 -74.87 140.02 1.00 7873
## time_cen:groupDM0.limbs -4.84 43.50 -90.02 80.63 1.00 10850
## time_cen:groupDM1.flaps -75.12 35.31 -145.53 -6.21 1.00 9525
## time_cen:groupDM1.limbs -103.14 43.41 -187.87 -16.43 1.00 10594
## Tail_ESS
## Intercept 6135
## time_cen 9155
## groupDM0.limbs 7652
## groupDM1.flaps 7317
## groupDM1.limbs 6860
## time_cen:groupDM0.limbs 9049
## time_cen:groupDM1.flaps 9589
## time_cen:groupDM1.limbs 8845
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 128.46 10.81 108.73 151.40 1.00 6797 8668
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(healing_model_int_rslope, effects='time_cen:group')
Open code
conditional_effects(healing_model_int_rslope, effects='group')
Again, the model converged well and shows sufficient ESS. Findings are the same but estimated uncertainty is larger, as expected.
2.3.2.5 Models diagnostic
Data density plot: Simulate data from models and compare it with actual data distributions
Open code
pla<-pp_check(healing_model_int, type='dens_overlay_grouped', ndraws = 50, group='group')
plb<-pp_check(healing_model_int_robust, type='dens_overlay_grouped', ndraws = 50, group='group')
plc<-pp_check(healing_model_int_rslope, type='dens_overlay_grouped', ndraws = 50, group='group')
plot_grid(pla, plb, plc, labels= c("A", "B", "C"), ncol = 1, nrow = 3)
The posterior predictive checks indicate good models’ fits, as the simulated data from the models closely resemble the observed data.
Scatter plot: are errors of the models larger for higher prediction?
Open code
pla<-pp_check(healing_model_int, type='scatter_avg')
plb<-pp_check(healing_model_int_robust, type='scatter_avg')
plc<-pp_check(healing_model_int_rslope, type='scatter_avg')
plot_grid(pla, plb, plc, labels= c("A", "B", "C"), ncol = 1, nrow = 3)
No sign of larger error for higher predictions
Min-max plot
Open code
pla<-pp_check(healing_model_int, type="stat_2d", stat = c("max", "min"), ndraws=200)
plb<-pp_check(healing_model_int_robust, type="stat_2d", stat = c("max", "min"), ndraws=200)
plc<-pp_check(healing_model_int_rslope, type="stat_2d", stat = c("max", "min"), ndraws=200)
plot_grid(pla, plb, plc, labels= c("A", "B", "C"), ncol = 1, nrow = 3)
While the models generally replicate the data well, a notable discrepancy arises with the prediction of lower minimum values compared to the actual dataset. Given that wound size is inherently non-negative, both Gaussian and robust regression models implicitly assume the possibility of negative outcomes, which does not align with our data characteristics. Despite this, the observed minimum values are not significantly higher than those predicted by the models, suggesting that the approximation remains reasonably accurate.
SD - mean plot
Open code
pla<-pp_check(healing_model_int, type="stat_2d", stat = c("mean", "sd"), ndraws=200)
plb<-pp_check(healing_model_int_robust, type="stat_2d", stat = c("mean", "sd"), ndraws=200)
plc<-pp_check(healing_model_int_rslope, type="stat_2d", stat = c("mean", "sd"), ndraws=200)
plot_grid(pla, plb, plc, labels= c("A", "B", "C"), ncol = 1, nrow = 3)
Look very well for all the models.
loo plot - calibartion
Open code
pla<-pp_check(healing_model_int, type = "loo_pit")
## Using all posterior draws for ppc type 'loo_pit' by default.
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
## Warning: '.fun' is deprecated.
## Use 'ppc_loo_pit_qq or ppc_loo_pit_overlay' instead.
## See help("Deprecated")
plb<-pp_check(healing_model_int_robust, type = "loo_pit")
## Using all posterior draws for ppc type 'loo_pit' by default.
## Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
## Warning: '.fun' is deprecated.
## Use 'ppc_loo_pit_qq or ppc_loo_pit_overlay' instead.
## See help("Deprecated")
plc<-pp_check(healing_model_int_rslope, type = "loo_pit")
## Using all posterior draws for ppc type 'loo_pit' by default.
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
## Warning: '.fun' is deprecated.
## Use 'ppc_loo_pit_qq or ppc_loo_pit_overlay' instead.
## See help("Deprecated")
plot_grid(pla, plb, plc, labels= c("A", "B", "C"), ncol = 1, nrow = 3)
Calibrations look well for all models.
loo plot - pareto values
Open code
par(mfrow=c(3,1))
plot(loo(healing_model_int))
plot(loo(healing_model_int_robust))
plot(loo(healing_model_int_rslope))
Robust model shows the best properties in terms of overly influential values.
2.3.2.6 Models comparison
There are some data points with Pareto >0.7, except for robust model. Let to provide moment_match and compare the models
Open code
healing_model_int <- add_criterion(healing_model_int,
criterion = "loo",
moment_match = TRUE)
healing_model_int_robust <- add_criterion(healing_model_int_robust,
criterion = "loo",
moment_match = TRUE)
healing_model_int_rslope <- add_criterion(healing_model_int_rslope,
criterion = "loo",
moment_match = TRUE)
loo_compare(healing_model_int, healing_model_int_robust, healing_model_int_rslope)
## elpd_diff se_diff
## healing_model_int_rslope 0.0 0.0
## healing_model_int_robust -6.3 3.8
## healing_model_int -7.5 4.7Although a random slope model offers some improvement in prediction accuracy, the enhancement is modest. Given the complexity of the random slope model and our dataset’s limited size (three observations per group), a simpler robust regression with a random intercept approach is preferred. This strategy acknowledges the potential underestimation of uncertainty, especially in group-specific slopes, but balances model complexity against our dataset’s constraints.
2.3.3 Vessel index model
As the vessel index also shows right-tailed distribution, generalized linear model with gamma distribution and log-link will be primarily used.
2.3.3.1 Model fit
Open code
vessel_model_int_gamma <- run_model(
brm(vessel_index ~ time_cen + group + time_cen:group + (1|Subject),
data = vessel_long,
prior = prior1,
chains = 2, iter = 8000, warmup = 2000, seed = 17,
control = list(adapt_delta = 0.99),
cores = 2,
family = Gamma(link = 'log') ),
'gitignore/output/brm/vessel_model_int_gamma', reuse = TRUE)
summary(vessel_model_int_gamma)
## Warning: There were 7 divergent transitions after warmup. Increasing
## adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Family: gamma
## Links: mu = log; shape = identity
## Formula: vessel_index ~ time_cen + group + time_cen:group + (1 | Subject)
## Data: vessel_long (Number of observations: 34)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Subject (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.17 0.12 0.01 0.48 1.00 2075 3496
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 4.16 0.14 3.86 4.43 1.00 5518
## time_cen -0.12 0.06 -0.24 -0.00 1.00 7167
## groupDM0.limbs 0.12 0.21 -0.28 0.55 1.00 5905
## groupDM1.flaps 0.10 0.19 -0.27 0.50 1.00 6019
## groupDM1.limbs 0.38 0.21 -0.01 0.81 1.00 5759
## time_cen:groupDM0.limbs 0.22 0.13 -0.04 0.50 1.00 7551
## time_cen:groupDM1.flaps 0.07 0.12 -0.17 0.32 1.00 8062
## time_cen:groupDM1.limbs -0.10 0.13 -0.36 0.18 1.00 8539
## Tail_ESS
## Intercept 4844
## time_cen 7844
## groupDM0.limbs 4984
## groupDM1.flaps 5035
## groupDM1.limbs 5048
## time_cen:groupDM0.limbs 6160
## time_cen:groupDM1.flaps 6236
## time_cen:groupDM1.limbs 6556
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 21.31 6.39 10.83 35.51 1.00 6021 7461
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(vessel_model_int_gamma, effects='time_cen:group')
Open code
conditional_effects(vessel_model_int_gamma, effects='group')
The model demonstrated good convergence and sufficient effective sample sizes (ESS), indicating reliable estimates. However, the vessel index variable presents challenges, primarily due to the lack of longitudinal data for most animals, resulting in significantly increased uncertainty in the estimates. Notably, the largest discrepancies are observed in the initial measurement. This considerable uncertainty makes it impossible to draw any conclusion from these results.
2.3.3.2 Posterior predictive check
Open code
pp_check(vessel_model_int_gamma, type='dens_overlay_grouped', ndraws = 50,group='group')
Open code
pp_check(vessel_model_int_gamma, type='scatter_avg')
Open code
pp_check(vessel_model_int_gamma, type="stat_2d", stat = c("max", "min"), ndraws=200)
Open code
pp_check(vessel_model_int_gamma, type="stat_2d", stat = c("mean", "sd"), ndraws=200)
Open code
pp_check(vessel_model_int_gamma, type = "loo_pit")
Open code
plot(loo(vessel_model_int_gamma))
Open code
loo(vessel_model_int_gamma)
##
## Computed from 12000 by 34 log-likelihood matrix
##
## Estimate SE
## elpd_loo -149.5 3.8
## p_loo 10.9 1.7
## looic 299.0 7.6
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 22 64.7% 1215
## (0.5, 0.7] (ok) 10 29.4% 873
## (0.7, 1] (bad) 2 5.9% 188
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.As above, PPC seems fine but with some overly influential values.
2.4 Results presentation
Model-based estimates and their uncertainty will be presented in the form of two figures and two tables per each of selected models (4 models are presented in total - one per outcome and additional random slope model for the wound area).
Figures
A) Visualizes the posterior distributions for differences between groups across three distinct measures:
- At the first measurement
- At the final measurement
- Based on the rate of change, highlighting how rapidly outcomes change in one group relative to another.
The posterior distributions will be presented to compare the following groups, arranged in the specified sequence (from top to bottom):
- Diabetes vs. non-DM in limbs
- Diabetes vs. non-DM in flaps
- Diabetes vs. non-DM, averaged over both limbs and flaps
- Limbs vs. flaps in non-diabetic subgroup
- Limbs vs. flaps in diabetic subgroup
- Limbs vs. flaps, averaged over both non-diabetic and diabetic
B) Visualizes the model-estimated change of the outcome, separately for each of 4 groups:
- Limb wounds in non-diabetic pigs (
Non-DM, limbs) - Flap wounds in non-diabetic pigs (
Non-DM, flaps) - Limb wounds in diabetic pigs (
Diabetic, limbs) - Flap wounds in diabetic pigs (
Diabetic, flaps)
Tables
Tables show estimates, 95% credible intervals (95% CI) and the probability of direction (PD) which is an index (0.5 to 1) representing the certainty that the effect goes in a particular direction, as described in Makowski et al., 2019
A) shows prediction of outcome separately for each of the 4 groups (see above) in terms of the following estimates: - (i) at the 1st measurement - (ii) at the final measurement - (iii) rate of change in in the outcome per a week
B) shows the estimated between-groups differences, for the same estimands as above
2.4.1 TcPO2
Posterior samples extraction
Open code
## original parameters (groups in middle time)
post_fix_tcpo <- as.data.frame(tcpo_model_int,
variable = c("b_Intercept","b_time_cen",
"b_groupDM0.limbs",
"b_groupDM1.flaps",
"b_groupDM1.limbs",
"b_time_cen:groupDM0.limbs",
"b_time_cen:groupDM1.flaps",
"b_time_cen:groupDM1.limbs"))
post_fix_tcpo <- post_fix_tcpo %>% mutate(
### posterior distributions of outcome across groups at middle time
DM0.flaps = b_Intercept,
DM0.limbs = b_Intercept + b_groupDM0.limbs,
DM1.flaps = b_Intercept + b_groupDM1.flaps,
DM1.limbs = b_Intercept + b_groupDM1.limbs,
### posterior distributions of slopes across groups
DM0.flaps_slope = b_time_cen,
DM0.limbs_slope = b_time_cen + post_fix_tcpo$'b_time_cen:groupDM0.limbs',
DM1.flaps_slope = b_time_cen + post_fix_tcpo$'b_time_cen:groupDM1.flaps',
DM1.limbs_slope = b_time_cen + post_fix_tcpo$'b_time_cen:groupDM1.limbs'
) %>% mutate(
### posterior distributions of outcome across group at the starting time
DM0.flaps_T0 = DM0.flaps - 2*DM0.flaps_slope,
DM0.limbs_T0 = DM0.limbs - 2*DM0.limbs_slope,
DM1.flaps_T0 = DM1.flaps - 2*DM1.flaps_slope,
DM1.limbs_T0 = DM1.limbs - 2*DM1.limbs_slope,
DM0.flaps_TF = DM0.flaps + 2*DM0.flaps_slope,
DM0.limbs_TF = DM0.limbs + 2*DM0.limbs_slope,
DM1.flaps_TF = DM1.flaps + 2*DM1.flaps_slope,
DM1.limbs_TF = DM1.limbs + 2*DM1.limbs_slope
) %>% mutate(
### difference: diabetes 0 vs. 1 in flaps
flaps.DMdif_TM = -DM0.flaps + DM1.flaps,
flaps.DMdif_T0 = -DM0.flaps_T0 + DM1.flaps_T0,
flaps.DMdif_TF = -DM0.flaps_TF + DM1.flaps_TF,
flaps.DMdif_slope = -DM0.flaps_slope + DM1.flaps_slope,
### difference: diabetes 0 vs. 1 in limbs
limbs.DMdif_TM = -DM0.limbs +DM1.limbs,
limbs.DMdif_T0 = -DM0.limbs_T0 +DM1.limbs_T0,
limbs.DMdif_TF = -DM0.limbs_TF +DM1.limbs_TF,
limbs.DMdif_slope = -DM0.limbs_slope +DM1.limbs_slope,
### difference in flaps vs. limbs in diabetes 0
DM0.sitedif_TM = -DM0.flaps + DM0.limbs,
DM0.sitedif_T0 = -DM0.flaps_T0 + DM0.limbs_T0,
DM0.sitedif_TF = -DM0.flaps_TF + DM0.limbs_TF,
DM0.sitedif_slope = -DM0.flaps_slope + DM0.limbs_slope,
### difference in flaps vs. limbs in diabetes 1
DM1.sitedif_TM = -DM1.flaps + DM1.limbs,
DM1.sitedif_T0 = -DM1.flaps_T0 + DM1.limbs_T0,
DM1.sitedif_TF = -DM1.flaps_TF + DM1.limbs_TF,
DM1.sitedif_slope = -DM1.flaps_slope + DM1.limbs_slope)
## difference: diabetes 0 vs. 1 averaged over both limbs and flaps
post_fix_tcpo$DMdif_TM <- rowMeans(
post_fix_tcpo[,c('flaps.DMdif_TM', 'limbs.DMdif_TM')]
)
post_fix_tcpo$DMdif_T0 <- rowMeans(
post_fix_tcpo[,c('flaps.DMdif_T0', 'limbs.DMdif_T0')]
)
post_fix_tcpo$DMdif_TF <- rowMeans(
post_fix_tcpo[,c('flaps.DMdif_TF', 'limbs.DMdif_TF')]
)
post_fix_tcpo$DMdif_slope <- rowMeans(
post_fix_tcpo[,c('flaps.DMdif_slope', 'limbs.DMdif_slope')]
)
## difference: flaps vs. limbs, averaged over both DM statuses
post_fix_tcpo$sitedif_TM <- rowMeans(
post_fix_tcpo[,c('DM0.sitedif_TM', 'DM1.sitedif_TM')]
)
post_fix_tcpo$sitedif_T0 <- rowMeans(
post_fix_tcpo[,c('DM0.sitedif_T0', 'DM1.sitedif_T0')]
)
post_fix_tcpo$sitedif_TF <- rowMeans(
post_fix_tcpo[,c('DM0.sitedif_TF', 'DM1.sitedif_TF')]
)
post_fix_tcpo$sitedif_slope <- rowMeans(
post_fix_tcpo[,c('DM0.sitedif_slope', 'DM1.sitedif_slope')]
)2.4.1.1 Posterior distributions of between-groups difference, SupplFig2
Following plot show posterior distribution for difference between the groups in terms of starting value, value at the end of the experiment and difference between predicted 1-week changes of the value. As the estimation was done with Gamma model with log link, the effects are shown in log-ratios, i.e. zero indicates null effect.
Open code
## general setting
suppl_fig_01 <- function(){
m <- matrix(c(1, 3, 4, 5,
2, 3, 4, 5), nrow = 2, ncol = 4, byrow = TRUE)
layout(mat = m,
heights = c(0.5, 0.5),
widths = c(0.04,rep(0.96/3,3)))
par(mgp=c(1.6,0.62,0))
par(mar=c(0,0,0,0))
plot(NULL, axes=FALSE,xlab="",ylab="",xlim=c(-1,1),ylim=c(-0.85,0.85))
text(0, 0 ,"Limbs vs flaps",
cex=1.3, font=3, col='blue',
xpd=TRUE, srt=90 )
lines(
c( .9 , .9),
c( -.7, .7),
lty = 2, col='blue' )
plot(NULL, axes=FALSE,xlab="",ylab="",xlim=c(-1,1),ylim=c(-0.85,0.85))
text(0, 0 ,"Diabetic vs non-diabetic", col='red3',
cex=1.3, font=3,
xpd=TRUE, srt=90)
lines(
c( .9 , .9),
c( -.55, .85),
lty = 2, col='red3' )
par(mar=c(2, 0, 1, 0))
par(mgp = c(0.7, 0.5, 0))
## T0 -------------------------------------------------------------------------------
xsc <- seq(-2, 3, 1)
ysc <- c(0, 6)
xscal <- max(xsc) - min(xsc)
plot(
NULL, axes=FALSE,
xlab="log-ratio at start",
ylab="",
cex.lab = 1.15,
xlim = c( min(xsc), max(xsc) ) ,
ylim = c( min(ysc), max(ysc) )
)
rect( min(xsc), 0, max(xsc) , max(ysc), col='grey91', border=NA)
for(i in 1:length(xsc) ) {
lines(c(xsc[i], xsc[i]), c(0, -0.02*max(ysc) ), xpd=TRUE)
lines(c(xsc[i], xsc[i]), c(0, ysc[2]), xpd=TRUE, col = 'white')
}
for(i in 1:length(xsc) ) {
text(
xsc[i], -0.04*max(ysc), xsc[i], xpd=TRUE
)
}
lines(c(0,0), c(0, max(ysc)), lty=2)
lines(c(min(xsc), max(xsc) ), c(0, 0), xpd=TRUE)
y <- 0.05*max(ysc)
ci(post_fix_tcpo, 'DMdif_T0', y, 'red3', xsc);y = y + 0.04*max(ysc)
cur(post_fix_tcpo, 'DMdif_T0', y, rgb(1, 0, 0, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix_tcpo, 'flaps.DMdif_T0', y, 'red3', xsc);y = y + 0.06*max(ysc)
ci(post_fix_tcpo, 'limbs.DMdif_T0', y, 'red3', xsc);y = y + 0.16*max(ysc)
ci(post_fix_tcpo, 'sitedif_T0', y, 'blue', xsc);y = y + 0.04*max(ysc)
cur(post_fix_tcpo, 'sitedif_T0', y, rgb(.1, 0.1, 1, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix_tcpo, 'DM1.sitedif_T0', y, 'blue', xsc);y = y + 0.06*max(ysc)
ci(post_fix_tcpo, 'DM0.sitedif_T0', y, 'blue', xsc)
text(min(xsc) + 0.2*xscal,
max(ysc)*0.975,
"non-DM subgroup", cex = 1.1)
text(min(xsc) + 0.2*xscal,
max(ysc)*0.915,
"diabetic subgroup", cex = 1.1)
text(min(xsc) + 0.17*xscal,
max(ysc)*0.66,
"average effect", cex = 1.1)
text(min(xsc) + 0.18*xscal,
max(ysc)*0.45,
"limbs subgroup", cex = 1.1)
text(min(xsc) + 0.18*xscal,
max(ysc)*0.38,
"flaps subgroup", cex = 1.1)
text(min(xsc) + 0.17*xscal,
max(ysc)*0.13,
"average effect", cex = 1.1)
text(min(xsc) + 0.8*xscal,
max(ysc)*0.7,
bquote("PD= "* .(
round(
p_dir(post_fix_tcpo$sitedif_T0, 'max', 0) , 2
)
)
),
col= 'blue' )
text(min(xsc) + 0.8*xscal,
max(ysc)*0.15,
bquote("PD= "* .(
round(
p_dir(post_fix_tcpo$DMdif_T0, 'max', 0) , 2
)
)
),
col= 'red3' )
## TF --------------------------------------------------------------------------------
xsc = seq(-2, 3, 1)
ysc = c(0, 6)
xscal <- max(xsc) - min(xsc)
plot(
NULL, axes=FALSE,
xlab="log-ratio at the end",
ylab="",
cex.lab = 1.15,
xlim = c( min(xsc), max(xsc) ) ,
ylim = c( min(ysc), max(ysc) )
)
rect( min(xsc), 0, max(xsc), max(ysc), col='grey91', border=NA)
for(i in 1:length(xsc) ) {
lines(c(xsc[i], xsc[i]), c(0, -0.02*max(ysc) ), xpd=TRUE)
lines(c(xsc[i], xsc[i]), c(0, ysc[2]), xpd=TRUE, col = 'white')
}
for(i in 1:length(xsc) ) {
text(
xsc[i], -0.04*max(ysc), xsc[i], xpd=TRUE
)
}
lines(c(0,0), c(0, max(ysc)), lty=2)
lines(c(min(xsc), max(xsc) ), c(0, 0), xpd=TRUE)
y <- 0.05*max(ysc)
ci(post_fix_tcpo, 'DMdif_TF', y, 'red3', xsc);y = y + 0.04*max(ysc)
cur(post_fix_tcpo, 'DMdif_TF', y, rgb(1, 0, 0, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix_tcpo, 'flaps.DMdif_TF', y, 'red3', xsc);y = y + 0.06*max(ysc)
ci(post_fix_tcpo, 'limbs.DMdif_TF', y, 'red3', xsc);y = y + 0.16*max(ysc)
ci(post_fix_tcpo, 'sitedif_TF', y, 'blue', xsc);y = y + 0.04*max(ysc)
cur(post_fix_tcpo, 'sitedif_TF', y, rgb(.1, 0.1, 1, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix_tcpo, 'DM1.sitedif_TF', y, 'blue', xsc);y = y + 0.06*max(ysc)
ci(post_fix_tcpo, 'DM0.sitedif_TF', y, 'blue', xsc)
text(min(xsc) + 0.8*xscal,
max(ysc)*0.7,
bquote("PD= "* .(
round(
p_dir(post_fix_tcpo$sitedif_TF, 'max', 0) , 3
)
)
),
col= 'blue' )
text(min(xsc) + 0.8*xscal,
max(ysc)*0.15,
bquote("PD= "* .(
round(
p_dir(post_fix_tcpo$DMdif_TF, 'max', 0) , 2
)
)
),
col= 'red3' )
## slope ----------------------------------------------------------------------------
xsc = round(seq(-0.8, 1, by = 0.4), 5)
ysc = c(0, 20)
xscal <- max(xsc) - min(xsc)
plot(
NULL, axes=FALSE,
xlab="log-ratio in magnitudes of change over a week",
ylab="",
cex.lab = 1.15,
xlim = c( min(xsc), max(xsc) ) ,
ylim = c( min(ysc), max(ysc) )
)
rect( min(xsc), 0, max(xsc) , max(ysc), col='grey91', border=NA)
for(i in 1:length(xsc) ) {
lines(c(xsc[i], xsc[i]), c(0, -0.02*max(ysc) ), xpd=TRUE)
lines(c(xsc[i], xsc[i]), c(0, ysc[2]), xpd=TRUE, col = 'white')
}
for(i in 1:length(xsc) ) {
text(
xsc[i], -0.04*max(ysc), xsc[i], xpd=TRUE
)
}
lines(c(0,0), c(0, max(ysc)), lty=2)
lines(c(min(xsc), max(xsc) ), c(0, 0), xpd=TRUE)
y <- 0.05*max(ysc)
ci(post_fix_tcpo, 'DMdif_slope', y, 'red3', xsc);y = y + 0.04*max(ysc)
cur(post_fix_tcpo, 'DMdif_slope', y, rgb(1, 0, 0, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix_tcpo, 'flaps.DMdif_slope', y, 'red3', xsc);y = y + 0.06*max(ysc)
ci(post_fix_tcpo, 'limbs.DMdif_slope', y, 'red3', xsc);y = y + 0.16*max(ysc)
ci(post_fix_tcpo, 'sitedif_slope', y, 'blue', xsc);y = y + 0.04*max(ysc)
cur(post_fix_tcpo, 'sitedif_slope', y, rgb(.1, 0.1, 1, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix_tcpo, 'DM1.sitedif_slope', y, 'blue', xsc);y = y + 0.06*max(ysc)
ci(post_fix_tcpo, 'DM0.sitedif_slope', y, 'blue', xsc)
text(min(xsc) + 0.8*xscal,
max(ysc)*0.7,
bquote("PD= "* .(
round(
p_dir(post_fix_tcpo$sitedif_slope, 'max', 0) , 2
)
)
),
col= 'blue' )
text(min(xsc) + 0.8*xscal,
max(ysc)*0.15,
bquote("PD= "* .(
round(
p_dir(post_fix_tcpo$DMdif_slope, 'max', 0) , 3
)
)
),
col= 'red3' )
}
suppl_fig_01()
The figure suggest that TcPO4 may restore faster in limbs, compared to flaps.
2.4.1.2 Model-estimated TcPO2 over time, Figure 1
Open code
time <- seq(0, 28, length=100)
time_z <- seq(-2, 2, length=100)
cip_flaps_DM0 <- data.frame()
for(i in 1:length(time_z)){
cip_flaps_DM0[1:12000,i] <- post_fix_tcpo$DM0.flaps + post_fix_tcpo$DM0.flaps_slope*time_z[i]}
cip_flaps_DM0 <- sapply(cip_flaps_DM0 , function(p) quantile(p, probs = c(0.025,0.975,0.5)))
cip_flaps_DM1 <- data.frame()
for(i in 1:length(time_z)){
cip_flaps_DM1[1:12000,i] <- post_fix_tcpo$DM1.flaps + post_fix_tcpo$DM1.flaps_slope*time_z[i]}
cip_flaps_DM1 <- sapply(cip_flaps_DM1 , function(p) quantile(p, probs = c(0.025,0.975,0.5)))
cip_limbs_DM0 <- data.frame()
for(i in 1:length(time_z)){
cip_limbs_DM0[1:12000,i] <- post_fix_tcpo$DM0.limbs + post_fix_tcpo$DM0.limbs_slope*time_z[i]}
cip_limbs_DM0 <- sapply(cip_limbs_DM0 , function(p) quantile(p, probs = c(0.025,0.975,0.5)))
cip_limbs_DM1 <- data.frame()
for(i in 1:length(time_z)){
cip_limbs_DM1[1:12000,i] <- post_fix_tcpo$DM1.limbs + post_fix_tcpo$DM1.limbs_slope*time_z[i]}
cip_limbs_DM1 <- sapply(cip_limbs_DM1 , function(p) quantile(p, probs = c(0.025,0.975,0.5)))
predict <- data.frame(
prediction =
unlist(c(
exp(cip_flaps_DM0[3,]), exp(cip_limbs_DM0[3,]), exp(cip_flaps_DM1[3,]), exp(cip_limbs_DM1[3,])
)),
cil =
unlist(c(
exp(cip_flaps_DM0[1,]), exp(cip_limbs_DM0[1,]), exp(cip_flaps_DM1[1,]), exp(cip_limbs_DM1[1,])
)),
ciu =
unlist(c(
exp(cip_flaps_DM0[2,]), exp(cip_limbs_DM0[2,]), exp(cip_flaps_DM1[2,]), exp(cip_limbs_DM1[2,])
)),
time = rep(time, 4), group = c(
rep("flaps_DM0",100), rep("limbs_DM0",100), rep("flaps_DM1",100), rep("limbs_DM1",100)
))
predict$group <- factor(predict$group, levels = c(
"flaps_DM0","limbs_DM0","flaps_DM1","limbs_DM1") )
cole <- c("#01AF40", "#0000FF","#FF0000",
"#FF00FF")
range <- c(0, 100)
tcPO_long <- tcPO_long %>%
mutate(group = fct_recode(group,
'flaps_DM0' = 'DM0.flaps',
'limbs_DM0' = 'DM0.limbs',
'flaps_DM1' = 'DM1.flaps',
'limbs_DM1' = 'DM1.limbs'
))
fig_01c <- ggplot() +
geom_line(data = tcPO_long,
aes(x = time_days, y = TcPO2, group = Site_nest),
alpha = 0.8, color = 'grey70') +
geom_line(data = predict, aes(x = time,
y = prediction,
group = group,
color = group), size=1.05) +
geom_ribbon(data = predict,
aes(x = time, ymin = cil, ymax = ciu, fill = group),
alpha = 0.25, color = NA) +
scale_color_manual(values = cole) +
scale_fill_manual(values = cole) +
labs(x = "Days of experiment", y = 'TcPO2 (mmHg)') +
facet_wrap(~group, ncol = 4, labeller = labeller(group = c("flaps_DM0" = "Non-DM, flaps",
"limbs_DM0" = "Non-DM, limbs",
"flaps_DM1" = "Diabetic, flaps",
"limbs_DM1" = "Diabetic, limbs"))) +
scale_x_continuous(breaks = c(0, 7, 14, 28))
fig_01c
2.4.1.3 TcPO2-relevant tables
Data preparation for the start of an experiment (T0)
Open code
quant <- post_fix_tcpo[, c(
'DM0.sitedif_T0',
'DM1.sitedif_T0',
'sitedif_T0',
'flaps.DMdif_T0',
'limbs.DMdif_T0',
'DMdif_T0' ) ]
qe <- data.frame()
for(i in 1:6){
qe[i,1:3] <- quantile(exp(quant[, i]),
probs = c(0.5, 0.025, 0.975)
)
}
qe <- round(qe, 2)
pd <- c()
for(i in 1:6) {
pd[i] <- p_dir(quant[,i], 'max', 0)
}
contrast = c('Limbs vs Flaps (non-diabetic subgroup)',
'Limbs vs Flaps (diabetic subgroup)',
'Limbs vs Flaps (average effect)',
'Diabetic vs non-DM (flaps subgroup)',
'Diabetic vs non-DM (limbs subgroup)',
'Diabetic vs non-DM (average effect)'
)
tcpo_T0_eff <- data.frame(
Ratio = c(qe[,1]),
CI_L = c(qe[,2]),
CI_U = c(qe[,3]),
PD = round(pd, 3)
)For the end of experiment (TF)
Open code
quant <- post_fix_tcpo[, c(
'DM0.sitedif_TF',
'DM1.sitedif_TF',
'sitedif_TF',
'flaps.DMdif_TF',
'limbs.DMdif_TF',
'DMdif_TF' ) ]
qe <- data.frame()
for(i in 1:6){
qe[i,1:3] <- quantile(exp(quant[, i]),
probs = c(0.5, 0.025, 0.975)
)
}
qe <- round(qe, 2)
pd <- c()
for(i in 1:6) {
pd[i] <- p_dir(quant[,i], 'max', 0)
}
tcpo_TF_eff <- data.frame(
#contrast,
Ratio = c(qe[,1]),
CI_L = c(qe[,2]),
CI_U = c(qe[,3]),
PD = round(pd, 3)
)Differences in slope of changes (fold-difference in 1-week fold-changes)
Open code
quant <- post_fix_tcpo[, c(
'DM0.sitedif_slope',
'DM1.sitedif_slope',
'sitedif_slope',
'flaps.DMdif_slope',
'limbs.DMdif_slope',
'DMdif_slope' ) ]
qe <- data.frame()
for(i in 1:6){
qe[i,1:3] <- quantile(exp(quant[, i]),
probs = c(0.5, 0.025, 0.975)
)
}
qe <- round(qe, 2)
pd <- c()
for(i in 1:6) {
pd[i] <- p_dir(quant[,i], 'max', 0)
}
tcpo_slope_eff <- data.frame(
#contrast,
Ratio = c(qe[,1]),
CI_L = c(qe[,2]),
CI_U = c(qe[,3]),
PD = round(pd, 3)
)
tot_eff <- cbind(tcpo_T0_eff, tcpo_TF_eff, tcpo_slope_eff)
row.names(tot_eff) <- contrast Estimating TcPO2 across groups
Open code
quant <- post_fix_tcpo[, c(
'DM0.flaps_T0',
'DM0.limbs_T0',
'DM1.flaps_T0',
'DM1.limbs_T0',
'limbs.DMdif_slope' ) ]
qe <- data.frame()
for(i in 1:4){
qe[i,1:3] <- quantile(exp(quant[, i]),
probs = c(0.5, 0.025, 0.975)
) }
qe <- round(qe, 2)
tcpo_T0_groupEst <- data.frame(
Estimate = c(qe[,1]),
CI_L = c(qe[,2]),
CI_U = c(qe[,3]) )
quant <- post_fix_tcpo[, c(
'DM0.flaps_TF',
'DM0.limbs_TF',
'DM1.flaps_TF',
'DM1.limbs_TF',
'limbs.DMdif_slope' ) ]
qe <- data.frame()
for(i in 1:4){
qe[i,1:3] <- quantile(exp(quant[, i]),
probs = c(0.5, 0.025, 0.975)
) }
qe <- round(qe, 2)
tcpo_TF_groupEst <- data.frame(
Estimate = c(qe[,1]),
CI_L = c(qe[,2]),
CI_U = c(qe[,3]) )
quant <- post_fix_tcpo[, c(
'DM0.flaps_slope',
'DM0.limbs_slope',
'DM1.flaps_slope',
'DM1.limbs_slope',
'limbs.DMdif_slope' ) ]
qe <- data.frame()
for(i in 1:4){
qe[i,1:3] <- quantile(exp(quant[, i]),
probs = c(0.5, 0.025, 0.975)
) }
qe <- round(qe, 2)
tcpo_slope_groupEst <- data.frame(
Slope = c(qe[,1]),
CI_L = c(qe[,2]),
CI_U = c(qe[,3]) )
tot_est <- cbind(tcpo_T0_groupEst, tcpo_TF_groupEst, tcpo_slope_groupEst)
row.names(tot_est) <- c("Non-DM flaps",
"Non-DM limbs",
"Diabetic flaps",
"Diabetic limbs") Printing the tables
Open code
table_01 <- kbl(tot_eff, caption =
"Table 1. Estimated differences in TcPO2 between groups at the start (A) and at the end (B) of the experiment, and difference between 1-week fold-change (C, slope). CI_L and CI_U: bounds of the 95% credible interval. PD: probability of direction. Estimates are based on a Bayesian hierarchical generalized linear model with a Gamma distribution and a log-link function.") %>%
kable_styling("striped",full_width = T) %>%
add_header_above(c(" " = 1, "(A) at the start" = 4,
"(B) at the end" = 4, "(C) difference in slope" = 4)) %>%
column_spec(1, width_min = '2.4in')
table_01| Ratio | CI_L | CI_U | PD | Ratio | CI_L | CI_U | PD | Ratio | CI_L | CI_U | PD | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Limbs vs Flaps (non-diabetic subgroup) | 2.96 | 1.08 | 8.21 | 0.982 | 2.40 | 0.72 | 7.39 | 0.924 | 0.94 | 0.68 | 1.31 | 0.631 |
| Limbs vs Flaps (diabetic subgroup) | 0.69 | 0.25 | 1.93 | 0.770 | 2.79 | 0.91 | 8.89 | 0.964 | 1.42 | 1.04 | 1.92 | 0.987 |
| Limbs vs Flaps (average effect) | 1.44 | 0.69 | 2.95 | 0.847 | 2.60 | 1.12 | 5.70 | 0.986 | 1.16 | 0.93 | 1.45 | 0.902 |
| Diabetic vs non-DM (flaps subgroup) | 3.14 | 1.22 | 8.02 | 0.989 | 1.10 | 0.35 | 2.95 | 0.574 | 0.76 | 0.58 | 1.01 | 0.973 |
| Diabetic vs non-DM (limbs subgroup) | 0.74 | 0.24 | 2.22 | 0.715 | 1.27 | 0.37 | 4.40 | 0.653 | 1.15 | 0.79 | 1.65 | 0.774 |
| Diabetic vs non-DM (average effect) | 1.52 | 0.72 | 3.18 | 0.875 | 1.19 | 0.51 | 2.59 | 0.663 | 0.94 | 0.75 | 1.17 | 0.715 |
Open code
suppl_table_01 <- kbl(tot_est, caption =
"Suppl. Table 1. Predicted TcPO2 levels at the start (A) and at the end (B) of the experiment, and the 1-week fold-change (C) across groups. CI_L and CI_U represent the bounds of the 95% credible interval. Estimates are derived from a Bayesian hierarchical generalized linear model with a Gamma distribution and a log-link function.") %>%
kable_styling("striped",full_width = T) %>%
add_header_above(c(" " = 1, "(A) TcPO2 at the start" = 3,
"(B) TcPO2 at the end" = 3, "(C) 1-week fold-change" = 3)) %>%
column_spec(1, width_min = '1.5in')
suppl_table_01| Estimate | CI_L | CI_U | Estimate | CI_L | CI_U | Slope | CI_L | CI_U | |
|---|---|---|---|---|---|---|---|---|---|
| Non-DM flaps | 3.97 | 2.05 | 7.81 | 12.79 | 6.42 | 30.42 | 1.35 | 1.10 | 1.67 |
| Non-DM limbs | 11.74 | 5.41 | 26.61 | 30.89 | 13.08 | 76.89 | 1.27 | 0.98 | 1.66 |
| Diabetic flaps | 12.53 | 6.21 | 24.59 | 14.13 | 6.82 | 29.34 | 1.03 | 0.86 | 1.23 |
| Diabetic limbs | 8.63 | 3.99 | 19.33 | 38.90 | 16.83 | 97.73 | 1.46 | 1.14 | 1.89 |
2.4.2 Wound area
Posterior draws extraction
Open code
## original parameters (groups in middle time)
post_fix_healing <- as.data.frame(healing_model_int_robust,
variable = c("b_Intercept","b_time_cen",
"b_groupDM0.limbs",
"b_groupDM1.flaps",
"b_groupDM1.limbs",
"b_time_cen:groupDM0.limbs",
"b_time_cen:groupDM1.flaps",
"b_time_cen:groupDM1.limbs"))
post_fix_healing <- post_fix_healing %>% mutate(
### posterior distributions of outcome across groups at middle time
DM0.flaps = b_Intercept,
DM0.limbs = b_Intercept + b_groupDM0.limbs,
DM1.flaps = b_Intercept + b_groupDM1.flaps,
DM1.limbs = b_Intercept + b_groupDM1.limbs,
### posterior distributions of slopes across groups
DM0.flaps_slope = b_time_cen,
DM0.limbs_slope = b_time_cen + post_fix_healing$'b_time_cen:groupDM0.limbs',
DM1.flaps_slope = b_time_cen + post_fix_healing$'b_time_cen:groupDM1.flaps',
DM1.limbs_slope = b_time_cen + post_fix_healing$'b_time_cen:groupDM1.limbs'
) %>% mutate(
### posterior distributions of outcome across group at the starting time
DM0.flaps_T0 = DM0.flaps - 1*DM0.flaps_slope,
DM0.limbs_T0 = DM0.limbs - 1*DM0.limbs_slope,
DM1.flaps_T0 = DM1.flaps - 1*DM1.flaps_slope,
DM1.limbs_T0 = DM1.limbs - 1*DM1.limbs_slope,
DM0.flaps_TF = DM0.flaps + 1*DM0.flaps_slope,
DM0.limbs_TF = DM0.limbs + 1*DM0.limbs_slope,
DM1.flaps_TF = DM1.flaps + 1*DM1.flaps_slope,
DM1.limbs_TF = DM1.limbs + 1*DM1.limbs_slope
) %>% mutate(
### difference: diabetes 0 vs. 1 in flaps
flaps.DMdif_TM = -DM0.flaps + DM1.flaps,
flaps.DMdif_T0 = -DM0.flaps_T0 + DM1.flaps_T0,
flaps.DMdif_TF = -DM0.flaps_TF + DM1.flaps_TF,
flaps.DMdif_slope = -DM0.flaps_slope + DM1.flaps_slope,
### difference: diabetes 0 vs. 1 in limbs
limbs.DMdif_TM = -DM0.limbs +DM1.limbs,
limbs.DMdif_T0 = -DM0.limbs_T0 +DM1.limbs_T0,
limbs.DMdif_TF = -DM0.limbs_TF +DM1.limbs_TF,
limbs.DMdif_slope = -DM0.limbs_slope +DM1.limbs_slope,
### difference in flaps vs. limbs in diabetes 0
DM0.sitedif_TM = -DM0.flaps + DM0.limbs,
DM0.sitedif_T0 = -DM0.flaps_T0 + DM0.limbs_T0,
DM0.sitedif_TF = -DM0.flaps_TF + DM0.limbs_TF,
DM0.sitedif_slope = -DM0.flaps_slope + DM0.limbs_slope,
### difference in flaps vs. limbs in diabetes 1
DM1.sitedif_TM = -DM1.flaps + DM1.limbs,
DM1.sitedif_T0 = -DM1.flaps_T0 + DM1.limbs_T0,
DM1.sitedif_TF = -DM1.flaps_TF + DM1.limbs_TF,
DM1.sitedif_slope = -DM1.flaps_slope + DM1.limbs_slope)
## difference: diabetes 0 vs. 1 averaged over both limbs and flaps
post_fix_healing$DMdif_TM <- rowMeans(
post_fix_healing[,c('flaps.DMdif_TM', 'limbs.DMdif_TM')]
)
post_fix_healing$DMdif_T0 <- rowMeans(
post_fix_healing[,c('flaps.DMdif_T0', 'limbs.DMdif_T0')]
)
post_fix_healing$DMdif_TF <- rowMeans(
post_fix_healing[,c('flaps.DMdif_TF', 'limbs.DMdif_TF')]
)
post_fix_healing$DMdif_slope <- rowMeans(
post_fix_healing[,c('flaps.DMdif_slope', 'limbs.DMdif_slope')]
)
## difference: flaps vs. limbs, averaged over both DM statuses
post_fix_healing$sitedif_TM <- rowMeans(
post_fix_healing[,c('DM0.sitedif_TM', 'DM1.sitedif_TM')]
)
post_fix_healing$sitedif_T0 <- rowMeans(
post_fix_healing[,c('DM0.sitedif_T0', 'DM1.sitedif_T0')]
)
post_fix_healing$sitedif_TF <- rowMeans(
post_fix_healing[,c('DM0.sitedif_TF', 'DM1.sitedif_TF')]
)
post_fix_healing$sitedif_slope <- rowMeans(
post_fix_healing[,c('DM0.sitedif_slope', 'DM1.sitedif_slope')]
)2.4.2.1 Posterior distributions of between-groups difference, SupplFig3
Open code
## general setting
suppl_fig_02 <- function(){
m <- matrix(c(1, 3, 4, 5,
2, 3, 4, 5), nrow = 2, ncol = 4, byrow = TRUE)
layout(mat = m,
heights = c(0.5, 0.5),
widths = c(0.04,rep(0.96/3,3)))
par(mgp=c(1.6,0.62,0))
par(mar=c(0,0,0,0))
plot(NULL, axes=FALSE,xlab="",ylab="",xlim=c(-1,1),ylim=c(-0.85,0.85))
text(0, 0 ,"Limbs vs flaps",
cex=1.3, font=3, col='blue',
xpd=TRUE, srt=90 )
lines(
c( .9 , .9),
c( -.7, .7),
lty = 2, col='blue' )
plot(NULL, axes=FALSE,xlab="",ylab="",xlim=c(-1,1),ylim=c(-0.85,0.85))
text(0, 0 ,"Diabetic vs non-DM", col='red3',
cex=1.3, font=3,
xpd=TRUE, srt=90)
lines(
c( .9 , .9),
c( -.55, .85),
lty = 2, col='red3' )
par(mar=c(2, 0, 1, 0))
par(mgp = c(0.7, 0.5, 0))
## T0 -------------------------------------------------------------------------------
xsc = seq(-600, 650, 200)
ysc = c(0, 0.04)
xscal <- max(xsc) - min(xsc)
plot(
NULL, axes=FALSE,
xlab="Difference in 1st measurement",
ylab="",
cex.lab = 1.15,
xlim = c( min(xsc), max(xsc) ) ,
ylim = c( min(ysc), max(ysc) )
)
rect( min(xsc), 0, max(xsc) , max(ysc), col='grey91', border=NA)
for(i in 1:length(xsc) ) {
lines(c(xsc[i], xsc[i]), c(0, -0.02*max(ysc) ), xpd=TRUE)
lines(c(xsc[i], xsc[i]), c(0, ysc[2]), xpd=TRUE, col = 'white')
}
for(i in 1:length(xsc) ) {
text(
xsc[i], -0.04*max(ysc), xsc[i], xpd=TRUE
)
}
lines(c(0,0), c(0, max(ysc)), lty=2)
lines(c(min(xsc), max(xsc) ), c(0, 0), xpd=TRUE)
y <- 0.05*max(ysc)
ci(post_fix_healing, 'DMdif_T0', y, 'red3', xsc);y = y + 0.04*max(ysc)
cur(post_fix_healing, 'DMdif_T0', y, rgb(1, 0, 0, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix_healing, 'flaps.DMdif_T0', y, 'red3', xsc);y = y + 0.06*max(ysc)
ci(post_fix_healing, 'limbs.DMdif_T0', y, 'red3', xsc);y = y + 0.16*max(ysc)
ci(post_fix_healing, 'sitedif_T0', y, 'blue', xsc);y = y + 0.04*max(ysc)
cur(post_fix_healing, 'sitedif_T0', y, rgb(.1, 0.1, 1, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix_healing, 'DM1.sitedif_T0', y, 'blue', xsc);y = y + 0.06*max(ysc)
ci(post_fix_healing, 'DM0.sitedif_T0', y, 'blue', xsc)
text(min(xsc) + 0.19*xscal,
max(ysc)*0.975,
"non-DM subgroup", cex = 1.1)
text(min(xsc) + 0.19*xscal,
max(ysc)*0.915,
"diabetic subgroup", cex = 1.1)
text(min(xsc) + 0.17*xscal,
max(ysc)*0.66,
"average effect", cex = 1.1)
text(min(xsc) + 0.17*xscal,
max(ysc)*0.45,
"limbs subgroup", cex = 1.1)
text(min(xsc) + 0.17*xscal,
max(ysc)*0.38,
"flaps subgroup", cex = 1.1)
text(min(xsc) + 0.17*xscal,
max(ysc)*0.13,
"average effect", cex = 1.1)
text(min(xsc) + 0.8*xscal,
max(ysc)*0.7,
bquote("PD= "* .(
round(
p_dir(post_fix_healing$sitedif_T0, 'max', 0) , 2
)
)
),
col= 'blue' )
text(min(xsc) + 0.8*xscal,
max(ysc)*0.15,
bquote("PD= "* .(
round(
p_dir(post_fix_healing$DMdif_T0, 'max', 0) , 2
)
)
),
col= 'red3' )
## TF --------------------------------------------------------------------------------
xsc = seq(-400, 650, 200)
ysc = c(0, .05)
xscal <- max(xsc) - min(xsc)
plot(
NULL, axes=FALSE,
xlab="Difference after a week",
ylab="",
cex.lab = 1.15,
xlim = c( min(xsc), max(xsc) ) ,
ylim = c( min(ysc), max(ysc) )
)
rect( min(xsc), 0, max(xsc) , max(ysc), col='grey91', border=NA)
for(i in 1:length(xsc) ) {
lines(c(xsc[i], xsc[i]), c(0, -0.02*max(ysc) ), xpd=TRUE)
lines(c(xsc[i], xsc[i]), c(0, ysc[2]), xpd=TRUE, col = 'white')
}
for(i in 1:length(xsc) ) {
text(
xsc[i], -0.04*max(ysc), xsc[i], xpd=TRUE
)
}
lines(c(0,0), c(0, max(ysc)), lty=2)
lines(c(min(xsc), max(xsc) ), c(0, 0), xpd=TRUE)
y <- 0.05*max(ysc)
ci(post_fix_healing, 'DMdif_TF', y, 'red3', xsc);y = y + 0.04*max(ysc)
cur(post_fix_healing, 'DMdif_TF', y, rgb(1, 0, 0, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix_healing, 'flaps.DMdif_TF', y, 'red3', xsc);y = y + 0.06*max(ysc)
ci(post_fix_healing, 'limbs.DMdif_TF', y, 'red3', xsc);y = y + 0.16*max(ysc)
ci(post_fix_healing, 'sitedif_TF', y, 'blue', xsc);y = y + 0.04*max(ysc)
cur(post_fix_healing, 'sitedif_TF', y, rgb(.1, 0.1, 1, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix_healing, 'DM1.sitedif_TF', y, 'blue', xsc);y = y + 0.06*max(ysc)
ci(post_fix_healing, 'DM0.sitedif_TF', y, 'blue', xsc)
text(min(xsc) + 0.8*xscal,
max(ysc)*0.7,
bquote("PD= "* .(
round(
p_dir(post_fix_healing$sitedif_TF, 'max', 0) , 3
)
)
),
col= 'blue' )
text(min(xsc) + 0.8*xscal,
max(ysc)*0.15,
bquote("PD= "* .(
round(
p_dir(post_fix_healing$DMdif_TF, 'max', 0) , 2
)
)
),
col= 'red3' )
## slope ----------------------------------------------------------------------------
xsc = round(seq(-300, 220, by = 100), 5)
ysc = c(0, 0.07)
xscal <- max(xsc) - min(xsc)
plot(
NULL, axes=FALSE,
xlab="Difference in change over time",
ylab="",
cex.lab = 1.15,
xlim = c( min(xsc), max(xsc) ) ,
ylim = c( min(ysc), max(ysc) )
)
rect( min(xsc), 0, max(xsc) , max(ysc), col='grey91', border=NA)
for(i in 1:length(xsc) ) {
lines(c(xsc[i], xsc[i]), c(0, -0.02*max(ysc) ), xpd=TRUE)
lines(c(xsc[i], xsc[i]), c(0, ysc[2]), xpd=TRUE, col = 'white')
}
for(i in 1:length(xsc) ) {
text(
xsc[i], -0.04*max(ysc), xsc[i], xpd=TRUE
)
}
lines(c(0,0), c(0, max(ysc)), lty=2)
lines(c(min(xsc), max(xsc) ), c(0, 0), xpd=TRUE)
y <- 0.05*max(ysc)
ci(post_fix_healing, 'DMdif_slope', y, 'red3', xsc);y = y + 0.04*max(ysc)
cur(post_fix_healing, 'DMdif_slope', y, rgb(1, 0, 0, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix_healing, 'flaps.DMdif_slope', y, 'red3', xsc);y = y + 0.06*max(ysc)
ci(post_fix_healing, 'limbs.DMdif_slope', y, 'red3', xsc);y = y + 0.16*max(ysc)
ci(post_fix_healing, 'sitedif_slope', y, 'blue', xsc);y = y + 0.04*max(ysc)
cur(post_fix_healing, 'sitedif_slope', y, rgb(.1, 0.1, 1, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix_healing, 'DM1.sitedif_slope', y, 'blue', xsc);y = y + 0.06*max(ysc)
ci(post_fix_healing, 'DM0.sitedif_slope', y, 'blue', xsc)
text(min(xsc) + 0.8*xscal,
max(ysc)*0.7,
bquote("PD= "* .(
round(
p_dir(post_fix_healing$sitedif_slope, 'max', 0) , 2
)
)
),
col= 'blue' )
text(min(xsc) + 0.8*xscal,
max(ysc)*0.15,
bquote("PD= "* .(
round(
p_dir(post_fix_healing$DMdif_slope, 'max', 0) , 3
)
)
),
col= 'red3' )
}
## Printing
suppl_fig_02()
2.4.2.2 Model-estimated wound area over time, Figure 2
Open code
time <- seq(0, 14, length=100)
time_z <- seq(-1, 1, length=100)
cip_flaps_DM0 <- data.frame()
for(i in 1:length(time_z)){
cip_flaps_DM0[1:12000,i] <- post_fix_healing$DM0.flaps + post_fix_healing$DM0.flaps_slope*time_z[i]}
cip_flaps_DM0 <- sapply(cip_flaps_DM0 , function(p) quantile(p, probs = c(0.025,0.975,0.5)))
cip_flaps_DM1 <- data.frame()
for(i in 1:length(time_z)){
cip_flaps_DM1[1:12000,i] <- post_fix_healing$DM1.flaps + post_fix_healing$DM1.flaps_slope*time_z[i]}
cip_flaps_DM1 <- sapply(cip_flaps_DM1 , function(p) quantile(p, probs = c(0.025,0.975,0.5)))
cip_limbs_DM0 <- data.frame()
for(i in 1:length(time_z)){
cip_limbs_DM0[1:12000,i] <- post_fix_healing$DM0.limbs + post_fix_healing$DM0.limbs_slope*time_z[i]}
cip_limbs_DM0 <- sapply(cip_limbs_DM0 , function(p) quantile(p, probs = c(0.025,0.975,0.5)))
cip_limbs_DM1 <- data.frame()
for(i in 1:length(time_z)){
cip_limbs_DM1[1:12000,i] <- post_fix_healing$DM1.limbs + post_fix_healing$DM1.limbs_slope*time_z[i]}
cip_limbs_DM1 <- sapply(cip_limbs_DM1 , function(p) quantile(p, probs = c(0.025,0.975,0.5)))
predict <- data.frame(
prediction =
unlist(c(
(cip_flaps_DM0[3,]), (cip_limbs_DM0[3,]), (cip_flaps_DM1[3,]), (cip_limbs_DM1[3,])
)),
cil =
unlist(c(
(cip_flaps_DM0[1,]), (cip_limbs_DM0[1,]), (cip_flaps_DM1[1,]), (cip_limbs_DM1[1,])
)),
ciu =
unlist(c(
(cip_flaps_DM0[2,]), (cip_limbs_DM0[2,]), (cip_flaps_DM1[2,]), (cip_limbs_DM1[2,])
)),
time = rep(time, 4), group = c(
rep("flaps_DM0",100), rep("limbs_DM0",100), rep("flaps_DM1",100), rep("limbs_DM1",100)
))
predict$group <- factor(predict$group, levels = c(
"flaps_DM0","limbs_DM0","flaps_DM1","limbs_DM1") )
cole <- c("#01AF40", "#0000FF","#FF0000",
"#FF00FF")
healing_long <- healing_long %>%
mutate(group = fct_recode(group,
'flaps_DM0' = 'DM0.flaps',
'limbs_DM0' = 'DM0.limbs',
'flaps_DM1' = 'DM1.flaps',
'limbs_DM1' = 'DM1.limbs'
))
fig_02c <- ggplot() +
geom_line(data = healing_long,
aes(x = time_days, y = healing, group = Site_nest),
alpha = 0.8, color = 'grey70') +
geom_line(data = predict, aes(x = time,
y = prediction,
group = group,
color = group), size=1.05) +
geom_ribbon(data = predict,
aes(x = time, ymin = cil, ymax = ciu, fill = group),
alpha = 0.25, color = NA) +
scale_color_manual(values = cole) +
scale_fill_manual(values = cole) +
labs(x = "Days of experiment", y = expression(Wound~area~(mm^2))) +
facet_wrap(~group, ncol = 4, labeller = labeller(group = c("flaps_DM0" = "Non-DM, flaps",
"limbs_DM0" = "Non-DM, limbs",
"flaps_DM1" = "Diabetic, flaps",
"limbs_DM1" = "Diabetic, limbs"))) +
scale_x_continuous(breaks=c(0, 7, 14))
fig_02c
2.4.3 Wound area, random slope
Posterior draws extraction
Open code
## original parameters (groups in middle time)
post_fix_healing_rslope <- as.data.frame(healing_model_int_rslope,
variable = c("b_Intercept","b_time_cen",
"b_groupDM0.limbs",
"b_groupDM1.flaps",
"b_groupDM1.limbs",
"b_time_cen:groupDM0.limbs",
"b_time_cen:groupDM1.flaps",
"b_time_cen:groupDM1.limbs"))
post_fix_healing_rslope <- post_fix_healing_rslope %>% mutate(
### posterior distributions of outcome across groups at middle time
DM0.flaps = b_Intercept,
DM0.limbs = b_Intercept + b_groupDM0.limbs,
DM1.flaps = b_Intercept + b_groupDM1.flaps,
DM1.limbs = b_Intercept + b_groupDM1.limbs,
### posterior distributions of slopes across groups
DM0.flaps_slope = b_time_cen,
DM0.limbs_slope = b_time_cen + post_fix_healing_rslope$'b_time_cen:groupDM0.limbs',
DM1.flaps_slope = b_time_cen + post_fix_healing_rslope$'b_time_cen:groupDM1.flaps',
DM1.limbs_slope = b_time_cen + post_fix_healing_rslope$'b_time_cen:groupDM1.limbs'
) %>% mutate(
### posterior distributions of outcome across group at the starting time
DM0.flaps_T0 = DM0.flaps - 1*DM0.flaps_slope,
DM0.limbs_T0 = DM0.limbs - 1*DM0.limbs_slope,
DM1.flaps_T0 = DM1.flaps - 1*DM1.flaps_slope,
DM1.limbs_T0 = DM1.limbs - 1*DM1.limbs_slope,
DM0.flaps_TF = DM0.flaps + 1*DM0.flaps_slope,
DM0.limbs_TF = DM0.limbs + 1*DM0.limbs_slope,
DM1.flaps_TF = DM1.flaps + 1*DM1.flaps_slope,
DM1.limbs_TF = DM1.limbs + 1*DM1.limbs_slope
) %>% mutate(
### difference: diabetes 0 vs. 1 in flaps
flaps.DMdif_TM = -DM0.flaps + DM1.flaps,
flaps.DMdif_T0 = -DM0.flaps_T0 + DM1.flaps_T0,
flaps.DMdif_TF = -DM0.flaps_TF + DM1.flaps_TF,
flaps.DMdif_slope = -DM0.flaps_slope + DM1.flaps_slope,
### difference: diabetes 0 vs. 1 in limbs
limbs.DMdif_TM = -DM0.limbs +DM1.limbs,
limbs.DMdif_T0 = -DM0.limbs_T0 +DM1.limbs_T0,
limbs.DMdif_TF = -DM0.limbs_TF +DM1.limbs_TF,
limbs.DMdif_slope = -DM0.limbs_slope +DM1.limbs_slope,
### difference in flaps vs. limbs in diabetes 0
DM0.sitedif_TM = -DM0.flaps + DM0.limbs,
DM0.sitedif_T0 = -DM0.flaps_T0 + DM0.limbs_T0,
DM0.sitedif_TF = -DM0.flaps_TF + DM0.limbs_TF,
DM0.sitedif_slope = -DM0.flaps_slope + DM0.limbs_slope,
### difference in flaps vs. limbs in diabetes 1
DM1.sitedif_TM = -DM1.flaps + DM1.limbs,
DM1.sitedif_T0 = -DM1.flaps_T0 + DM1.limbs_T0,
DM1.sitedif_TF = -DM1.flaps_TF + DM1.limbs_TF,
DM1.sitedif_slope = -DM1.flaps_slope + DM1.limbs_slope)
## difference: diabetes 0 vs. 1 averaged over both limbs and flaps
post_fix_healing_rslope$DMdif_TM <- rowMeans(
post_fix_healing_rslope[,c('flaps.DMdif_TM', 'limbs.DMdif_TM')]
)
post_fix_healing_rslope$DMdif_T0 <- rowMeans(
post_fix_healing_rslope[,c('flaps.DMdif_T0', 'limbs.DMdif_T0')]
)
post_fix_healing_rslope$DMdif_TF <- rowMeans(
post_fix_healing_rslope[,c('flaps.DMdif_TF', 'limbs.DMdif_TF')]
)
post_fix_healing_rslope$DMdif_slope <- rowMeans(
post_fix_healing_rslope[,c('flaps.DMdif_slope', 'limbs.DMdif_slope')]
)
## difference: flaps vs. limbs, averaged over both DM statuses
post_fix_healing_rslope$sitedif_TM <- rowMeans(
post_fix_healing_rslope[,c('DM0.sitedif_TM', 'DM1.sitedif_TM')]
)
post_fix_healing_rslope$sitedif_T0 <- rowMeans(
post_fix_healing_rslope[,c('DM0.sitedif_T0', 'DM1.sitedif_T0')]
)
post_fix_healing_rslope$sitedif_TF <- rowMeans(
post_fix_healing_rslope[,c('DM0.sitedif_TF', 'DM1.sitedif_TF')]
)
post_fix_healing_rslope$sitedif_slope <- rowMeans(
post_fix_healing_rslope[,c('DM0.sitedif_slope', 'DM1.sitedif_slope')]
)Preparing data for tables
Start of an experiment (T0)
Open code
quant <- post_fix_healing_rslope[, c(
'DM0.sitedif_T0',
'DM1.sitedif_T0',
'sitedif_T0',
'flaps.DMdif_T0',
'limbs.DMdif_T0',
'DMdif_T0' ) ]
qe <- data.frame()
for(i in 1:6){
qe[i,1:3] <- quantile(quant[, i],
probs = c(0.5, 0.025, 0.975)
)
}
qe <- round(qe)
pd <- c()
for(i in 1:6) {
pd[i] <- p_dir(quant[,i], 'max', 0)
}
contrast = c('Limbs vs Flaps (non-diabetic subgroup)',
'Limbs vs Flaps (diabetic subgroup)',
'Limbs vs Flaps (average effect)',
'Diabetic vs non-DM (flaps subgroup)',
'Diabetic vs non-DM(limbs subgroup)',
'Diabetic vs non-DM (average effect)'
)
tcpo_T0_eff <- data.frame(
Diff = c(qe[,1]),
CI_L = c(qe[,2]),
CI_U = c(qe[,3]),
PD = round(pd, 3)
)At the end of experiment (TF)
Open code
quant <- post_fix_healing_rslope[, c(
'DM0.sitedif_TF',
'DM1.sitedif_TF',
'sitedif_TF',
'flaps.DMdif_TF',
'limbs.DMdif_TF',
'DMdif_TF' ) ]
qe <- data.frame()
for(i in 1:6){
qe[i,1:3] <- quantile(quant[, i],
probs = c(0.5, 0.025, 0.975)
)
}
qe <- round(qe)
pd <- c()
for(i in 1:6) {
pd[i] <- p_dir(quant[,i], 'max', 0)
}
tcpo_TF_eff <- data.frame(
Diff= c(qe[,1]),
CI_L = c(qe[,2]),
CI_U = c(qe[,3]),
PD = round(pd, 3)
)Differences in slope of changes (fold-difference in 1-week fold-changes)
Open code
quant <- post_fix_healing_rslope[, c(
'DM0.sitedif_slope',
'DM1.sitedif_slope',
'sitedif_slope',
'flaps.DMdif_slope',
'limbs.DMdif_slope',
'DMdif_slope' ) ]
qe <- data.frame()
for(i in 1:6){
qe[i,1:3] <- quantile(quant[, i],
probs = c(0.5, 0.025, 0.975)
)
}
qe <- round(qe)
pd <- c()
for(i in 1:6) {
pd[i] <- p_dir(quant[,i], 'max', 0)
}
tcpo_slope_eff <- data.frame(
Diff = c(qe[,1]),
CI_L = c(qe[,2]),
CI_U = c(qe[,3]),
PD = round(pd, 3)
)
tot_eff <- cbind(tcpo_T0_eff, tcpo_TF_eff, tcpo_slope_eff)
row.names(tot_eff) <- contrast Estimating wound area across groups
Open code
quant <- post_fix_healing_rslope[, c(
'DM0.flaps_T0',
'DM0.limbs_T0',
'DM1.flaps_T0',
'DM1.limbs_T0',
'limbs.DMdif_slope' ) ]
qe <- data.frame()
for(i in 1:4){
qe[i,1:3] <- quantile(quant[, i],
probs = c(0.5, 0.025, 0.975)
) }
qe <- round(qe)
tcpo_T0_groupEst <- data.frame(
Estimate = c(qe[,1]),
CI_L = c(qe[,2]),
CI_U = c(qe[,3]) )
quant <- post_fix_healing_rslope[, c(
'DM0.flaps_TF',
'DM0.limbs_TF',
'DM1.flaps_TF',
'DM1.limbs_TF',
'limbs.DMdif_slope' ) ]
qe <- data.frame()
for(i in 1:4){
qe[i,1:3] <- quantile(quant[, i],
probs = c(0.5, 0.025, 0.975)
) }
qe <- round(qe)
tcpo_TF_groupEst <- data.frame(
Estimate = c(qe[,1]),
CI_L = c(qe[,2]),
CI_U = c(qe[,3]) )
quant <- post_fix_healing_rslope[, c(
'DM0.flaps_slope',
'DM0.limbs_slope',
'DM1.flaps_slope',
'DM1.limbs_slope',
'limbs.DMdif_slope' ) ]
qe <- data.frame()
for(i in 1:4){
qe[i,1:3] <- quantile(quant[, i],
probs = c(0.5, 0.025, 0.975)
) }
qe <- round(qe)
tcpo_slope_groupEst <- data.frame(
Slope = c(qe[,1]),
CI_L = c(qe[,2]),
CI_U = c(qe[,3]) )
tot_est <- cbind(tcpo_T0_groupEst, tcpo_TF_groupEst, tcpo_slope_groupEst)
row.names(tot_est) <- c("Non-DM flaps",
"Non-DM limbs",
"Diabetic flaps",
"Diabetic limbs") 2.4.3.1 Tables
Printing the tables
Open code
suppl_table_03 <- kbl(tot_eff, caption =
"Suppl. Table 3. Estimated differences in wound area between groups on the start (A, D0) and at th end (B, D14) of experiment, and difference between the rate of change per a week (C, slope). CI_L and CI_U: bounds of 95% credible interval. PD: probability of direction. Estimates are based on Bayesian hierarchical linear model, modelling both random intercept and slope and their covariance") %>%
kable_styling("striped",full_width = T) %>%
add_header_above(c(" " = 1, "(A) on the start" = 4,
"(B) at the end" = 4, "(C) difference in slope" = 4)) %>%
column_spec(1, width_min = '2.2in')
suppl_table_04 <- kbl(tot_est, caption =
"Suppl. Table 4. Predicted wound area on the start (A, day 0) and at the end (B, 14th day) of experiment, and 1-week change, reflecting healing (C, change per week) across groups. CI_L and CI_U: bounds of 95% credible interval. Estimates are based on Bayesian hierarchical linear model, modelling both random intercept and slope and their covariance") %>%
kable_styling("striped",full_width = T) %>%
add_header_above(c(" " = 1, "(A) Wound area on the start" = 3,
"(B) Wound area at the end" = 3, "(C) 1-week change" = 3)) %>%
column_spec(1, width_min = '1.5in')
suppl_table_03| Diff | CI_L | CI_U | PD | Diff | CI_L | CI_U | PD | Diff | CI_L | CI_U | PD | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Limbs vs Flaps (non-diabetic subgroup) | -153 | -307 | 4 | 0.972 | -143 | -331 | 52 | 0.931 | 5 | -133 | 143 | 0.529 |
| Limbs vs Flaps (diabetic subgroup) | 131 | -34 | 297 | 0.943 | 76 | -129 | 286 | 0.783 | -28 | -177 | 129 | 0.653 |
| Limbs vs Flaps (average effect) | -10 | -125 | 104 | 0.574 | -33 | -174 | 109 | 0.691 | -12 | -111 | 93 | 0.598 |
| Diabetic vs non-DM (flaps subgroup) | -8 | -154 | 136 | 0.548 | -130 | -302 | 60 | 0.924 | -61 | -187 | 76 | 0.831 |
| Diabetic vs non-DM(limbs subgroup) | 276 | 95 | 448 | 0.998 | 94 | -128 | 314 | 0.810 | -90 | -251 | 70 | 0.884 |
| Diabetic vs non-DM (average effect) | 134 | 15 | 246 | 0.985 | -18 | -156 | 129 | 0.604 | -76 | -174 | 32 | 0.932 |
Open code
suppl_table_04| Estimate | CI_L | CI_U | Estimate | CI_L | CI_U | Slope | CI_L | CI_U | |
|---|---|---|---|---|---|---|---|---|---|
| Non-DM flaps | 489 | 391 | 594 | 335 | 200 | 458 | -77 | -175 | 14 |
| Non-DM limbs | 336 | 211 | 465 | 191 | 38 | 344 | -73 | -183 | 39 |
| Diabetic flaps | 481 | 373 | 588 | 206 | 72 | 346 | -138 | -238 | -33 |
| Diabetic limbs | 614 | 486 | 736 | 283 | 130 | 441 | -166 | -273 | -50 |
2.4.4 Vessel index
Posterior draws extraction
Open code
## original parameters (groups in middle time)
post_fix <- as.data.frame(vessel_model_int_gamma,
variable = c("b_Intercept","b_time_cen",
"b_groupDM0.limbs",
"b_groupDM1.flaps",
"b_groupDM1.limbs",
"b_time_cen:groupDM0.limbs",
"b_time_cen:groupDM1.flaps",
"b_time_cen:groupDM1.limbs"))
post_fix <- post_fix %>% mutate(
### posterior distributions of outcome across groups at middle time
DM0.flaps = b_Intercept,
DM0.limbs = b_Intercept + b_groupDM0.limbs,
DM1.flaps = b_Intercept + b_groupDM1.flaps,
DM1.limbs = b_Intercept + b_groupDM1.limbs,
### posterior distributions of slopes across groups
DM0.flaps_slope = b_time_cen,
DM0.limbs_slope = b_time_cen + post_fix$'b_time_cen:groupDM0.limbs',
DM1.flaps_slope = b_time_cen + post_fix$'b_time_cen:groupDM1.flaps',
DM1.limbs_slope = b_time_cen + post_fix$'b_time_cen:groupDM1.limbs'
) %>% mutate(
### posterior distributions of outcome across group at the starting time
DM0.flaps_T0 = DM0.flaps - 1.5*DM0.flaps_slope,
DM0.limbs_T0 = DM0.limbs - 1.5*DM0.limbs_slope,
DM1.flaps_T0 = DM1.flaps - 1.5*DM1.flaps_slope,
DM1.limbs_T0 = DM1.limbs - 1.5*DM1.limbs_slope,
DM0.flaps_TF = DM0.flaps + 1.5*DM0.flaps_slope,
DM0.limbs_TF = DM0.limbs + 1.5*DM0.limbs_slope,
DM1.flaps_TF = DM1.flaps + 1.5*DM1.flaps_slope,
DM1.limbs_TF = DM1.limbs + 1.5*DM1.limbs_slope
) %>% mutate(
### difference: diabetes 0 vs. 1 in flaps
flaps.DMdif_TM = -DM0.flaps + DM1.flaps,
flaps.DMdif_T0 = -DM0.flaps_T0 + DM1.flaps_T0,
flaps.DMdif_TF = -DM0.flaps_TF + DM1.flaps_TF,
flaps.DMdif_slope = -DM0.flaps_slope + DM1.flaps_slope,
### difference: diabetes 0 vs. 1 in limbs
limbs.DMdif_TM = -DM0.limbs +DM1.limbs,
limbs.DMdif_T0 = -DM0.limbs_T0 +DM1.limbs_T0,
limbs.DMdif_TF = -DM0.limbs_TF +DM1.limbs_TF,
limbs.DMdif_slope = -DM0.limbs_slope +DM1.limbs_slope,
### difference in flaps vs. limbs in diabetes 0
DM0.sitedif_TM = -DM0.flaps + DM0.limbs,
DM0.sitedif_T0 = -DM0.flaps_T0 + DM0.limbs_T0,
DM0.sitedif_TF = -DM0.flaps_TF + DM0.limbs_TF,
DM0.sitedif_slope = -DM0.flaps_slope + DM0.limbs_slope,
### difference in flaps vs. limbs in diabetes 1
DM1.sitedif_TM = -DM1.flaps + DM1.limbs,
DM1.sitedif_T0 = -DM1.flaps_T0 + DM1.limbs_T0,
DM1.sitedif_TF = -DM1.flaps_TF + DM1.limbs_TF,
DM1.sitedif_slope = -DM1.flaps_slope + DM1.limbs_slope)
## difference: diabetes 0 vs. 1 averaged over both limbs and flaps
post_fix$DMdif_TM <- rowMeans(
post_fix[,c('flaps.DMdif_TM', 'limbs.DMdif_TM')]
)
post_fix$DMdif_T0 <- rowMeans(
post_fix[,c('flaps.DMdif_T0', 'limbs.DMdif_T0')]
)
post_fix$DMdif_TF <- rowMeans(
post_fix[,c('flaps.DMdif_TF', 'limbs.DMdif_TF')]
)
post_fix$DMdif_slope <- rowMeans(
post_fix[,c('flaps.DMdif_slope', 'limbs.DMdif_slope')]
)
## difference: flaps vs. limbs, averaged over both DM statuses
post_fix$sitedif_TM <- rowMeans(
post_fix[,c('DM0.sitedif_TM', 'DM1.sitedif_TM')]
)
post_fix$sitedif_T0 <- rowMeans(
post_fix[,c('DM0.sitedif_T0', 'DM1.sitedif_T0')]
)
post_fix$sitedif_TF <- rowMeans(
post_fix[,c('DM0.sitedif_TF', 'DM1.sitedif_TF')]
)
post_fix$sitedif_slope <- rowMeans(
post_fix[,c('DM0.sitedif_slope', 'DM1.sitedif_slope')]
)2.4.4.1 Posterior distributions of between-groups difference, SupplFig3
Open code
## general setting
suppl_fig_03 <- function(){
m <- matrix(c(1, 3, 4, 5,
2, 3, 4, 5), nrow = 2, ncol = 4, byrow = TRUE)
layout(mat = m,
heights = c(0.5, 0.5),
widths = c(0.04,rep(0.96/3,3)))
par(mgp=c(1.6,0.62,0))
par(mar=c(0,0,0,0))
plot(NULL, axes=FALSE,xlab="",ylab="",xlim=c(-1,1),ylim=c(-0.85,0.85))
text(0, 0 ,"Limbs vs flaps",
cex=1.3, font=3, col='blue',
xpd=TRUE, srt=90 )
lines(
c( .9 , .9),
c( -.7, .7),
lty = 2, col='blue' )
plot(NULL, axes=FALSE,xlab="",ylab="",xlim=c(-1,1),ylim=c(-0.85,0.85))
text(0, 0 ,"Diabetic vs non-diabetic", col='red3',
cex=1.3, font=3,
xpd=TRUE, srt=90)
lines(
c( .9 , .9),
c( -.55, .85),
lty = 2, col='red3' )
par(mar=c(2, 0, 1, 0))
par(mgp = c(0.7, 0.5, 0))
## T0 -------------------------------------------------------------------------------
xsc = seq(-2, 3, 1)
ysc = c(0, 10)
xscal <- max(xsc) - min(xsc)
plot(
NULL, axes=FALSE,
xlab="log-ratio at start",
ylab="",
cex.lab = 1.15,
xlim = c( min(xsc), max(xsc) ) ,
ylim = c( min(ysc), max(ysc) )
)
rect( min(xsc), 0, max(xsc) , max(ysc), col='grey91', border=NA)
for(i in 1:length(xsc) ) {
lines(c(xsc[i], xsc[i]), c(0, -0.02*max(ysc) ), xpd=TRUE)
lines(c(xsc[i], xsc[i]), c(0, ysc[2]), xpd=TRUE, col = 'white')
}
for(i in 1:length(xsc) ) {
text(
xsc[i], -0.04*max(ysc), xsc[i], xpd=TRUE
)
}
lines(c(0,0), c(0, max(ysc)), lty=2)
lines(c(min(xsc), max(xsc) ), c(0, 0), xpd=TRUE)
y <- 0.05*max(ysc)
ci(post_fix, 'DMdif_T0', y, 'red3', xsc);y = y + 0.04*max(ysc)
cur(post_fix, 'DMdif_T0', y, rgb(1, 0, 0, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix, 'flaps.DMdif_T0', y, 'red3', xsc);y = y + 0.06*max(ysc)
ci(post_fix, 'limbs.DMdif_T0', y, 'red3', xsc);y = y + 0.16*max(ysc)
ci(post_fix, 'sitedif_T0', y, 'blue', xsc);y = y + 0.04*max(ysc)
cur(post_fix, 'sitedif_T0', y, rgb(.1, 0.1, 1, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix, 'DM1.sitedif_T0', y, 'blue', xsc);y = y + 0.06*max(ysc)
ci(post_fix, 'DM0.sitedif_T0', y, 'blue', xsc)
text(min(xsc) + 0.2*xscal,
max(ysc)*0.975,
"non-DM subgroup", cex = 1.1)
text(min(xsc) + 0.2*xscal,
max(ysc)*0.915,
"diabetic subgroup", cex = 1.1)
text(min(xsc) + 0.17*xscal,
max(ysc)*0.66,
"average effect", cex = 1.1)
text(min(xsc) + 0.18*xscal,
max(ysc)*0.45,
"limbs subgroup", cex = 1.1)
text(min(xsc) + 0.18*xscal,
max(ysc)*0.38,
"flaps subgroup", cex = 1.1)
text(min(xsc) + 0.17*xscal,
max(ysc)*0.13,
"average effect", cex = 1.1)
text(min(xsc) + 0.8*xscal,
max(ysc)*0.7,
bquote("PD= "* .(
round(
p_dir(post_fix$sitedif_T0, 'max', 0) , 2
)
)
),
col= 'blue' )
text(min(xsc) + 0.8*xscal,
max(ysc)*0.15,
bquote("PD= "* .(
round(
p_dir(post_fix$DMdif_T0, 'max', 0) , 3
)
)
),
col= 'red3' )
## TM --------------------------------------------------------------------------------
xsc = seq(-2, 2, 1)
ysc = c(0, 10)
xscal <- max(xsc) - min(xsc)
plot(
NULL, axes=FALSE,
xlab="log-ratio at the end",
ylab="",
cex.lab = 1.15,
xlim = c( min(xsc), max(xsc) ) ,
ylim = c( min(ysc), max(ysc) )
)
rect( min(xsc), 0, max(xsc), max(ysc), col='grey91', border=NA)
for(i in 1:length(xsc) ) {
lines(c(xsc[i], xsc[i]), c(0, -0.02*max(ysc) ), xpd=TRUE)
lines(c(xsc[i], xsc[i]), c(0, ysc[2]), xpd=TRUE, col = 'white')
}
for(i in 1:length(xsc) ) {
text(
xsc[i], -0.04*max(ysc), xsc[i], xpd=TRUE
)
}
lines(c(0,0), c(0, max(ysc)), lty=2)
lines(c(min(xsc), max(xsc) ), c(0, 0), xpd=TRUE)
y <- 0.05*max(ysc)
ci(post_fix, 'DMdif_TF', y, 'red3', xsc);y = y + 0.04*max(ysc)
cur(post_fix, 'DMdif_TF', y, rgb(1, 0, 0, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix, 'flaps.DMdif_TF', y, 'red3', xsc);y = y + 0.06*max(ysc)
ci(post_fix, 'limbs.DMdif_TF', y, 'red3', xsc);y = y + 0.16*max(ysc)
ci(post_fix, 'sitedif_TF', y, 'blue', xsc);y = y + 0.04*max(ysc)
cur(post_fix, 'sitedif_TF', y, rgb(.1, 0.1, 1, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix, 'DM1.sitedif_TF', y, 'blue', xsc);y = y + 0.06*max(ysc)
ci(post_fix, 'DM0.sitedif_TF', y, 'blue', xsc)
text(min(xsc) + 0.8*xscal,
max(ysc)*0.7,
bquote("PD= "* .(
round(
p_dir(post_fix$sitedif_TF, 'max', 0) , 3
)
)
),
col= 'blue' )
text(min(xsc) + 0.8*xscal,
max(ysc)*0.15,
bquote("PD= "* .(
round(
p_dir(post_fix$DMdif_TF, 'max', 0) , 2
)
)
),
col= 'red3' )
## slope ----------------------------------------------------------------------------
xsc = round(seq(-1, 1, by = 0.5), 5)
ysc = c(0, 20)
xscal <- max(xsc) - min(xsc)
plot(
NULL, axes=FALSE,
xlab="log-ratio in magnitudes of change over a week",
ylab="",
cex.lab = 1.15,
xlim = c( min(xsc), max(xsc) ) ,
ylim = c( min(ysc), max(ysc) )
)
rect( min(xsc), 0, max(xsc) , max(ysc), col='grey91', border=NA)
for(i in 1:length(xsc) ) {
lines(c(xsc[i], xsc[i]), c(0, -0.02*max(ysc) ), xpd=TRUE)
lines(c(xsc[i], xsc[i]), c(0, ysc[2]), xpd=TRUE, col = 'white')
}
for(i in 1:length(xsc) ) {
text(
xsc[i], -0.04*max(ysc), xsc[i], xpd=TRUE
)
}
lines(c(0,0), c(0, max(ysc)), lty=2)
lines(c(min(xsc), max(xsc) ), c(0, 0), xpd=TRUE)
y <- 0.05*max(ysc)
ci(post_fix, 'DMdif_slope', y, 'red3', xsc);y = y + 0.04*max(ysc)
cur(post_fix, 'DMdif_slope', y, rgb(1, 0, 0, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix, 'flaps.DMdif_slope', y, 'red3', xsc);y = y + 0.06*max(ysc)
ci(post_fix, 'limbs.DMdif_slope', y, 'red3', xsc);y = y + 0.16*max(ysc)
ci(post_fix, 'sitedif_slope', y, 'blue', xsc);y = y + 0.04*max(ysc)
cur(post_fix, 'sitedif_slope', y, rgb(.1, 0.1, 1, alpha=0.4));y = y + 0.27*max(ysc)
ci(post_fix, 'DM1.sitedif_slope', y, 'blue', xsc);y = y + 0.06*max(ysc)
ci(post_fix, 'DM0.sitedif_slope', y, 'blue', xsc)
text(min(xsc) + 0.8*xscal,
max(ysc)*0.7,
bquote("PD= "* .(
round(
p_dir(post_fix$sitedif_slope, 'max', 0) , 2
)
)
),
col= 'blue' )
text(min(xsc) + 0.8*xscal,
max(ysc)*0.15,
bquote("PD= "* .(
round(
p_dir(post_fix$DMdif_slope, 'max', 0) , 3
)
)
),
col= 'red3' )
}
suppl_fig_03()
The figure does not suggest any reliable between-group difference in terms of final values of vessel index or rate of change.
2.4.4.3 Model-estimated vessel index over time, Figure 3
Open code
time <- seq(7, 28, length=100)
time_z <- seq(-1.5, 1.5, length=100)
cip_flaps_DM0 <- data.frame()
for(i in 1:length(time_z)){
cip_flaps_DM0[1:12000,i] <- post_fix$DM0.flaps + post_fix$DM0.flaps_slope*time_z[i]}
cip_flaps_DM0 <- sapply(cip_flaps_DM0 , function(p) quantile(p, probs = c(0.025,0.975,0.5)))
cip_flaps_DM1 <- data.frame()
for(i in 1:length(time_z)){
cip_flaps_DM1[1:12000,i] <- post_fix$DM1.flaps + post_fix$DM1.flaps_slope*time_z[i]}
cip_flaps_DM1 <- sapply(cip_flaps_DM1 , function(p) quantile(p, probs = c(0.025,0.975,0.5)))
cip_limbs_DM0 <- data.frame()
for(i in 1:length(time_z)){
cip_limbs_DM0[1:12000,i] <- post_fix$DM0.limbs + post_fix$DM0.limbs_slope*time_z[i]}
cip_limbs_DM0 <- sapply(cip_limbs_DM0 , function(p) quantile(p, probs = c(0.025,0.975,0.5)))
cip_limbs_DM1 <- data.frame()
for(i in 1:length(time_z)){
cip_limbs_DM1[1:12000,i] <- post_fix$DM1.limbs + post_fix$DM1.limbs_slope*time_z[i]}
cip_limbs_DM1 <- sapply(cip_limbs_DM1 , function(p) quantile(p, probs = c(0.025,0.975,0.5)))
predict <- data.frame(
prediction =
unlist(c(
exp(cip_flaps_DM0[3,]), exp(cip_limbs_DM0[3,]), exp(cip_flaps_DM1[3,]), exp(cip_limbs_DM1[3,])
)),
cil =
unlist(c(
exp(cip_flaps_DM0[1,]), exp(cip_limbs_DM0[1,]), exp(cip_flaps_DM1[1,]), exp(cip_limbs_DM1[1,])
)),
ciu =
unlist(c(
exp(cip_flaps_DM0[2,]), exp(cip_limbs_DM0[2,]), exp(cip_flaps_DM1[2,]), exp(cip_limbs_DM1[2,])
)),
time = rep(time, 4), group = c(
rep("flaps_DM0",100), rep("limbs_DM0",100),rep("flaps_DM1",100), rep("limbs_DM1",100)
))
predict$group <- factor(predict$group, levels = c(
"flaps_DM0","limbs_DM0","flaps_DM1","limbs_DM1") )
cole <- c("#01AF40", "#0000FF","#FF0000",
"#FF00FF")
vessel_long <- vessel_long %>%
mutate(group = fct_recode(group,
'flaps_DM0' = 'DM0.flaps',
'limbs_DM0' = 'DM0.limbs',
'flaps_DM1' = 'DM1.flaps',
'limbs_DM1' = 'DM1.limbs'
))
fig_03c <- ggplot() +
geom_point(data = vessel_long,
aes(x = time_days, y = vessel_index, group = Site_nest),
alpha = 0.8, color = 'grey70') +
geom_line(data = predict, aes(x = time,
y = prediction,
group = group,
color = group), size=1.05) +
geom_ribbon(data = predict,
aes(x = time, ymin = cil, ymax = ciu, fill = group),
alpha = 0.25, color = NA) +
scale_color_manual(values = cole) +
scale_fill_manual(values = cole) +
labs(x = "Days of experiment", y = 'Vessel index') +
facet_wrap(~group, ncol = 4,
labeller = labeller(group = c("flaps_DM0" = "Non-DM, flaps",
"limbs_DM0" = "Non-DM, limbs",
"flaps_DM1" = "Diabetic, flaps",
"limbs_DM1" = "Diabetic, limbs"))) +
scale_x_continuous(breaks=c(7, 14, 28)) +
guides(color = guide_legend(override.aes = list(size = 4.6, linewidth=1), nrow=1)) +
theme(legend.position= 'no')
# ,
# legend.key.size = unit(0.7, 'cm'),
# legend.title = element_text(size=11),
# legend.text = element_text(size=10))
fig_03c
3 Findings description
3.1 TcPO2
Figure 1A shows TcPO2 trajectories in individual wounds over time. Wounds in non-diabetic animals on flaps maintained low TcPO2 levels throughout the whole experiment, whereas TcPO2 was markedly increasing in limbs. Initially, TcPO2 levels in control wounds differed markedly from those in occluded wounds, but limb wounds gradually matched control levels, unlike flap wounds, which remained distinct through day 28 (Fig. 1B) particularly in non-diabetic animals.
The results from the Bayesian hierarchical generalized linear model corroborate these findings (Fig. 1B, Table 1, Suppl. Fig. 2, Suppl. Table 1). The model highlights a significant effect of wound location, demonstrating that limb wounds exhibit a 2.6-fold higher TcPO2 value compared to flap wounds on day 28. Specifically, for non-diabetic and diabetic subgroups, the final TcPO2 values were approximately 2.4-fold and 2.8-fold greater in limb wounds, respectively (Table 1). The model estimates a 98.6% probability that limb location results in enhanced TcPO2 levels by day 28. The model found no reliable evidence of diabetes affecting TcPO2 levels (Table 1, Suppl. Fig. 2).
In conclusion, limb wounds seems linked to higher TcPO2 levels a month after wound formation, indicating better healing process for limb wounds.
3.1.1 Figure 1
Open code
fig_01c <- fig_01c + theme(legend.position = "none")
plot_grid(fig_01b, fig_01c, ncol=1, rel_heights = c(2.7, 1), labels=c("A", "B"))
3.1.2 Table 1
Open code
table_01| Ratio | CI_L | CI_U | PD | Ratio | CI_L | CI_U | PD | Ratio | CI_L | CI_U | PD | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Limbs vs Flaps (non-diabetic subgroup) | 2.96 | 1.08 | 8.21 | 0.982 | 2.40 | 0.72 | 7.39 | 0.924 | 0.94 | 0.68 | 1.31 | 0.631 |
| Limbs vs Flaps (diabetic subgroup) | 0.69 | 0.25 | 1.93 | 0.770 | 2.79 | 0.91 | 8.89 | 0.964 | 1.42 | 1.04 | 1.92 | 0.987 |
| Limbs vs Flaps (average effect) | 1.44 | 0.69 | 2.95 | 0.847 | 2.60 | 1.12 | 5.70 | 0.986 | 1.16 | 0.93 | 1.45 | 0.902 |
| Diabetic vs non-DM (flaps subgroup) | 3.14 | 1.22 | 8.02 | 0.989 | 1.10 | 0.35 | 2.95 | 0.574 | 0.76 | 0.58 | 1.01 | 0.973 |
| Diabetic vs non-DM (limbs subgroup) | 0.74 | 0.24 | 2.22 | 0.715 | 1.27 | 0.37 | 4.40 | 0.653 | 1.15 | 0.79 | 1.65 | 0.774 |
| Diabetic vs non-DM (average effect) | 1.52 | 0.72 | 3.18 | 0.875 | 1.19 | 0.51 | 2.59 | 0.663 | 0.94 | 0.75 | 1.17 | 0.715 |
3.1.3 Supplementary figure 2
Open code
suppl_fig_01()
3.1.4 Supplementary table 1
Open code
suppl_table_01| Estimate | CI_L | CI_U | Estimate | CI_L | CI_U | Slope | CI_L | CI_U | |
|---|---|---|---|---|---|---|---|---|---|
| Non-DM flaps | 3.97 | 2.05 | 7.81 | 12.79 | 6.42 | 30.42 | 1.35 | 1.10 | 1.67 |
| Non-DM limbs | 11.74 | 5.41 | 26.61 | 30.89 | 13.08 | 76.89 | 1.27 | 0.98 | 1.66 |
| Diabetic flaps | 12.53 | 6.21 | 24.59 | 14.13 | 6.82 | 29.34 | 1.03 | 0.86 | 1.23 |
| Diabetic limbs | 8.63 | 3.99 | 19.33 | 38.90 | 16.83 | 97.73 | 1.46 | 1.14 | 1.89 |
3.2 Wound area
Wound area declined over time (D0 to D14 post-wound formation) in all groups. Control (unoccluded) wounds showed results comparable to the experimental wounds (Fig. 2A). Diabetic wounds exhibited surprisingly steeper declining trends, indicating higher healing velocity (Fig. 2A).
The Bayesian model supports this observation (Fig. 2B, Table 2, Suppl. Fig. 3). Healing velocity — defined as the weekly rate of wound area reduction — showed a notable difference between groups. The rate of change per week was -74 mm2 in non-diabetic pigs and -158 mm2 in diabetic pigs, leading to a difference of 84 mm2 per week in healing velocity rate between the diabetic and non-diabetic pigs (95% CI: 23 to 144; Suppl. Table 2). Within the subgroups, diabetic pigs exhibited a greater healing velocity by 69 mm2 per week in flap wounds (CI: 1 to 138) and 99 mm2 per week in limb wounds (-3 to 199) compared to their non-diabetic counterparts (Table 2, Suppl. Fig. 3). Remarkably, the model suggests a 99.7% probability that diabetes enhances healing velocity, i.e. fasten wound closure over time.
For limb wounds, the initially larger wound area in diabetic pigs complicates the interpretation of results (Fig. 2C, Table 2). However, flap wounds started with similar sizes across both groups, yet data suggest diabetes still promotes quicker healing in diabetic flap wounds, with a 97.6% model-estimated probability of faster wound closure (Table 2).
Although these results seem relatively robust, sensitivity analysis advises caution (Suppl. Tables 3 and 4). A complex model incorporating a random slope yielded a smaller model-based probability of 93.2% for diabetes facilitating faster healing (Suppl. Table 3), suggesting uncertainty.
3.2.0.1 Figure 2
Open code
fig_02c <- fig_02c + theme(legend.position = "none")
plot_grid(fig_02b, fig_02c, ncol=1, rel_heights = c(2.7, 1), labels=c("A", "B"))
3.2.0.2 Table 2
Open code
table_02| Diff | CI_L | CI_U | PD | Diff | CI_L | CI_U | PD | Diff | CI_L | CI_U | PD | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Limbs vs Flaps (non-diabetic subgroup) | -153 | -286 | -19 | 0.988 | -153 | -288 | -12 | 0.984 | 1 | -86 | 87 | 0.505 |
| Limbs vs Flaps (diabetic subgroup) | 134 | -6 | 269 | 0.970 | 74 | -67 | 215 | 0.856 | -29 | -116 | 60 | 0.745 |
| Limbs vs Flaps (average effect) | -10 | -107 | 83 | 0.590 | -39 | -137 | 57 | 0.789 | -15 | -75 | 47 | 0.685 |
| Diabetic vs non-DM (flaps subgroup) | -1 | -115 | 115 | 0.512 | -139 | -255 | -20 | 0.987 | -69 | -138 | -1 | 0.976 |
| Diabetic vs non-DM (limbs subgroup) | 283 | 128 | 437 | 1.000 | 88 | -76 | 242 | 0.863 | -99 | -199 | 3 | 0.971 |
| Diabetic vs non-DM (average effect) | 141 | 46 | 237 | 0.997 | -25 | -124 | 71 | 0.700 | -84 | -144 | -23 | 0.997 |
3.2.0.3 Supplementary figure 3
Open code
suppl_fig_02()
3.2.0.4 Supplementary table 2
Open code
suppl_table_02| Estimate | CI_L | CI_U | Estimate | CI_L | CI_U | Slope | CI_L | CI_U | |
|---|---|---|---|---|---|---|---|---|---|
| Non-DM flaps | 486 | 407 | 569 | 338 | 255 | 424 | -74 | -123 | -24 |
| Non-DM limbs | 334 | 223 | 441 | 185 | 77 | 298 | -74 | -147 | -1 |
| Diabetic flaps | 485 | 402 | 570 | 199 | 116 | 285 | -143 | -191 | -94 |
| Diabetic limbs | 617 | 508 | 725 | 273 | 160 | 386 | -172 | -244 | -99 |
3.2.0.5 Supplementary table 3
Open code
suppl_table_03| Diff | CI_L | CI_U | PD | Diff | CI_L | CI_U | PD | Diff | CI_L | CI_U | PD | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Limbs vs Flaps (non-diabetic subgroup) | -153 | -307 | 4 | 0.972 | -143 | -331 | 52 | 0.931 | 5 | -133 | 143 | 0.529 |
| Limbs vs Flaps (diabetic subgroup) | 131 | -34 | 297 | 0.943 | 76 | -129 | 286 | 0.783 | -28 | -177 | 129 | 0.653 |
| Limbs vs Flaps (average effect) | -10 | -125 | 104 | 0.574 | -33 | -174 | 109 | 0.691 | -12 | -111 | 93 | 0.598 |
| Diabetic vs non-DM (flaps subgroup) | -8 | -154 | 136 | 0.548 | -130 | -302 | 60 | 0.924 | -61 | -187 | 76 | 0.831 |
| Diabetic vs non-DM(limbs subgroup) | 276 | 95 | 448 | 0.998 | 94 | -128 | 314 | 0.810 | -90 | -251 | 70 | 0.884 |
| Diabetic vs non-DM (average effect) | 134 | 15 | 246 | 0.985 | -18 | -156 | 129 | 0.604 | -76 | -174 | 32 | 0.932 |
3.2.0.6 Supplementary table 4
Open code
suppl_table_04| Estimate | CI_L | CI_U | Estimate | CI_L | CI_U | Slope | CI_L | CI_U | |
|---|---|---|---|---|---|---|---|---|---|
| Non-DM flaps | 489 | 391 | 594 | 335 | 200 | 458 | -77 | -175 | 14 |
| Non-DM limbs | 336 | 211 | 465 | 191 | 38 | 344 | -73 | -183 | 39 |
| Diabetic flaps | 481 | 373 | 588 | 206 | 72 | 346 | -138 | -238 | -33 |
| Diabetic limbs | 614 | 486 | 736 | 283 | 130 | 441 | -166 | -273 | -50 |
3.3 Vessel index
Data suggest different trends in the vessel index over time, with an increasing trend for wounds on limbs of non-diabetic pigs and a decline in flaps. However, this may be influenced by differing initial values (Fig. 3). Given that the data are not longitudinal (measured once per wound) and initial values vary significantly across groups (Fig. 3), interpreting the results is challenging. Congruently, the Bayesian model shows a large uncertainty in estimates and does not indicate any reliable difference between groups (Table 3, Suppl. Table 5).
3.3.1 Figure 3
Open code
fig_03c
3.3.2 Table 3
Open code
table_03| Ratio | CI_L | CI_U | PD | Ratio | CI_L | CI_U | PD | Ratio | CI_L | CI_U | PD | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Limbs vs Flaps (non-DM subgroup) | 0.81 | 0.48 | 1.42 | 0.811 | 1.54 | 0.87 | 2.98 | 0.946 | 1.24 | 0.96 | 1.65 | 0.957 |
| Limbs vs Flaps (diabetic subgroup) | 1.69 | 0.95 | 3.01 | 0.966 | 1.04 | 0.54 | 2.03 | 0.554 | 0.85 | 0.62 | 1.17 | 0.879 |
| Limbs vs Flaps (average effect) | 1.17 | 0.79 | 1.74 | 0.816 | 1.27 | 0.82 | 2.03 | 0.887 | 1.03 | 0.84 | 1.28 | 0.621 |
| Diabetic vs non-DM (flaps subgroup) | 0.99 | 0.60 | 1.68 | 0.512 | 1.20 | 0.71 | 2.16 | 0.791 | 1.07 | 0.84 | 1.37 | 0.727 |
| Diabetic vs non-DM (limbs subgroup) | 2.08 | 1.13 | 3.84 | 0.986 | 0.81 | 0.40 | 1.68 | 0.758 | 0.73 | 0.52 | 1.03 | 0.968 |
| Diabetic vs non-DM (average effect) | 1.44 | 0.97 | 2.16 | 0.966 | 0.99 | 0.63 | 1.60 | 0.527 | 0.88 | 0.72 | 1.10 | 0.902 |
3.3.3 Suppl. Table 5
Open code
suppl_table_05| Estimate | CI_L | CI_U | Estimate | CI_L | CI_U | Slope | CI_L | CI_U | |
|---|---|---|---|---|---|---|---|---|---|
| Non-DM flaps | 76.93 | 54.94 | 107.24 | 53.89 | 37.57 | 73.29 | 0.89 | 0.79 | 1.00 |
| Non-DM limbs | 62.01 | 40.78 | 98.28 | 83.48 | 50.71 | 141.62 | 1.10 | 0.86 | 1.40 |
| Diabetic flaps | 76.51 | 51.31 | 114.67 | 64.87 | 41.26 | 102.93 | 0.95 | 0.76 | 1.18 |
| Diabetic limbs | 129.19 | 84.86 | 203.73 | 67.37 | 41.24 | 113.00 | 0.80 | 0.63 | 1.02 |
3.3.4 Conclusion
TcPO2 appears to be restored faster in limb wounds compared to flaps, while, surprisingly, healing velocity is better in diabetic animals, potentially due to insulin treatment and the short duration of diabetes. Analysis of the vessel index does not provide robust evidence for the effect of diabetes or wound location.
Generally, flap wounds in non-diabetic animals exhibited the worst endpoint values across the three outcomes, suggesting this group may serve as the most suitable model for chronic diabetic wounds. Inducing of diabetes does not appear to reinforce chronicity of the wound in our pig models.
Given the small sample size, caution in interpretation is warranted. However, our publicly available data invite future Bayesian updating and meta-analysis, contributing to a deeper understanding of wound healing mechanisms and potential interventions.
4 Software versions
Open code
sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=cs_CZ.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=cs_CZ.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=cs_CZ.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=cs_CZ.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Prague
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] kableExtra_1.4.0 ggbeeswarm_0.7.2 bayesplot_1.8.1 corrplot_0.92
## [5] arm_1.13-1 lme4_1.1-35.1 MASS_7.3-60 janitor_2.2.0
## [9] projpred_2.8.0 brms_2.16.3 Rcpp_1.0.12 glmnet_4.1-8
## [13] Matrix_1.6-5 boot_1.3-30 cowplot_1.1.1 pROC_1.18.0
## [17] mgcv_1.9-1 nlme_3.1-163 openxlsx_4.2.5 flextable_0.9.5
## [21] sjPlot_2.8.15 car_3.1-2 carData_3.0-5 gtsummary_1.7.2
## [25] emmeans_1.7.2 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
## [29] purrr_1.0.2 readr_2.1.2 tidyr_1.3.1 tibble_3.2.1
## [33] ggplot2_3.5.0 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] shinythemes_1.2.0 splines_4.3.3 later_1.3.0
## [4] cellranger_1.1.0 xts_0.12.1 reprex_2.0.1
## [7] lifecycle_1.0.4 StanHeaders_2.21.0-7 processx_3.8.3
## [10] lattice_0.22-5 insight_0.19.8 crosstalk_1.2.0
## [13] backports_1.4.1 magrittr_2.0.3 rmarkdown_2.26
## [16] yaml_2.3.5 httpuv_1.6.5 zip_2.2.0
## [19] askpass_1.1 pkgbuild_1.3.1 DBI_1.1.2
## [22] minqa_1.2.4 lubridate_1.8.0 multcomp_1.4-18
## [25] abind_1.4-5 rvest_1.0.2 TH.data_1.1-0
## [28] tensorA_0.36.2 sandwich_3.0-1 gdtools_0.3.7
## [31] inline_0.3.19 crul_1.4.0 performance_0.10.9
## [34] bridgesampling_1.1-2 svglite_2.1.0 codetools_0.2-19
## [37] DT_0.17 xml2_1.3.3 tidyselect_1.2.0
## [40] shape_1.4.6.1 farver_2.1.1 ggeffects_1.5.0
## [43] httpcode_0.3.0 matrixStats_1.2.0 stats4_4.3.3
## [46] base64enc_0.1-3 broom.helpers_1.14.0 jsonlite_1.8.8
## [49] ellipsis_0.3.2 ggridges_0.5.3 survival_3.5-8
## [52] iterators_1.0.14 systemfonts_1.0.4 foreach_1.5.2
## [55] tools_4.3.3 ragg_1.2.1 glue_1.7.0
## [58] gridExtra_2.3 xfun_0.42 distributional_0.4.0
## [61] loo_2.4.1 withr_3.0.0 fastmap_1.1.1
## [64] fansi_1.0.6 shinyjs_2.1.0 openssl_1.4.6
## [67] callr_3.7.5 digest_0.6.34 R6_2.5.1
## [70] mime_0.12 estimability_1.3 textshaping_0.3.6
## [73] colorspace_2.1-0 gtools_3.9.2 markdown_1.12
## [76] threejs_0.3.3 utf8_1.2.4 generics_0.1.2
## [79] fontLiberation_0.1.0 data.table_1.15.2 prettyunits_1.1.1
## [82] httr_1.4.2 htmlwidgets_1.6.4 pkgconfig_2.0.3
## [85] dygraphs_1.1.1.6 gtable_0.3.4 rsconnect_0.8.25
## [88] htmltools_0.5.7 fontBitstreamVera_0.1.1 scales_1.3.0
## [91] posterior_1.2.0 snakecase_0.11.1 knitr_1.45
## [94] rstudioapi_0.13 tzdb_0.4.0 reshape2_1.4.4
## [97] uuid_1.0-3 coda_0.19-4 checkmate_2.0.0
## [100] curl_4.3.2 nloptr_2.0.0 zoo_1.8-9
## [103] sjlabelled_1.2.0 vipor_0.4.7 parallel_4.3.3
## [106] miniUI_0.1.1.1 pillar_1.9.0 grid_4.3.3
## [109] vctrs_0.6.5 shinystan_2.5.0 promises_1.2.0.1
## [112] dbplyr_2.1.1 xtable_1.8-4 beeswarm_0.4.0
## [115] evaluate_0.23 mvtnorm_1.1-3 cli_3.6.2
## [118] compiler_4.3.3 rlang_1.1.3 crayon_1.5.2
## [121] rstantools_2.1.1 labeling_0.4.2 modelr_0.1.8
## [124] ps_1.7.6 plyr_1.8.9 sjmisc_2.8.9
## [127] fs_1.6.3 stringi_1.8.3 rstan_2.21.3
## [130] viridisLite_0.4.2 assertthat_0.2.1 munsell_0.5.0
## [133] colourpicker_1.1.1 Brobdingnag_1.2-7 bayestestR_0.13.2
## [136] fontquiver_0.2.1 sjstats_0.18.2 hms_1.1.1
## [139] gfonts_0.2.0 shiny_1.8.0 highr_0.9
## [142] haven_2.4.3 igraph_1.5.1 gt_0.10.1
## [145] broom_1.0.5 RcppParallel_5.1.7 readxl_1.3.1
## [148] officer_0.6.5References
Bürkner, Paul-Christian. 2017. “Brms: An r Package for Bayesian Multilevel Models Using Stan” 80. https://doi.org/10.18637/jss.v080.i01.
———. 2018. “Advanced Bayesian Multilevel Modeling with the R Package Brms.” The R Journal 10 (1): 395–411. https://doi.org/10.32614/RJ-2018-017.
Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. Chapman; Hall/CRC. https://doi.org/10.1201/b16018.
Makowski, Dominique, Mattan S. Ben-Shachar, S. H. Annabel Chen, and Daniel Lüdecke. 2019. “Indices of Effect Existence and Significance in the Bayesian Framework.” Frontiers in Psychology 10 (December). https://doi.org/10.3389/fpsyg.2019.02767.
McElreath, Richard. 2020. Statistical Rethinking. Chapman; Hall/CRC. https://doi.org/10.1201/9780429029608.
R Core Team. 2023. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Roy, Sashwati, Savita Khanna, Cameron Rink, Sabyasachi Biswas, and Chandan K. Sen. 2008. “Characterization of the Acute Temporal Changes in Excisional Murine Cutaneous Wound Inflammation by Screening of the Wound-Edge Transcriptome.” Physiological Genomics 34 (2): 162–84. https://doi.org/10.1152/physiolgenomics.00045.2008.
Stan Development Team. 2021. “RStan: The R Interface to Stan.” https://mc-stan.org/.
Vehtari, Aki, Jonah Gabry, Mans Magnusson, Yuling Yao, Paul-Christian Bürkner, Topi Paananen, and Andrew Gelman. 2020. “Loo: Efficient Leave-One-Out Cross-Validation and WAIC for Bayesian Models.” https://mc-stan.org/loo/.
Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2016. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27 (5): 1413–32. https://doi.org/10.1007/s11222-016-9696-4.